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PREFACE 

Despite the existence of numerous authoritative books in the field of orbital 
mechanics, the author has felt the need for a book which places emphasis on the con- 
ditions encountered with actively controlled satellites and spaceships rather than on 
observation and analysis of the passive heavenly bodies treated in classical astron- 
omy. Present-day space research, including the use of computers, has made much of 
the material in previous books outmoded; less emphasis is now placed on closed-form 
solutions^and more on iterative techniques. It is also apparent that a greater empha- 
sis on the basic formulas has become necessary. The problem of relative motion 
between two vehicles, which was rarely encountered in classical astronomy, has 
become a routine operational matter today and deserves consideration appropriate to 
its present importance. Likewise, the atmospheric entry problem for controlled 
vehicles is new and should be treated. The author has endeavored to use a direct 
approach to the problems at hand, where possible, and has sought to avoid the cum- 
bersome methods employed by some of the previous works which were dictated in 
their day by a need to avoid as much arithmetic manipulation as possible. This book 
is intended, however, not as a competitor but as a supplement to other books. 

Although the material covered here is not altogether new, it is scattered in the 
research literature and is not always readily accessible. The purpose of this book is 
to provide the reader with a reasonably complete primer on modern dynamic orbit 
theory. To this end many simple steps, which would seem obvious to the accom- 
plished professional in orbital mechanics or flight dynamics, have been retained. For 
people already familiar with the subject, it is hoped that the book will provide a use- 
ful reference handbook. It should provide a convenient guide for the solution of 
practical problems encountered in the analysis of space missions or the design of 
satellites and space vehicles. 

To those who are not so accomplished in the field it is hoped that the essentially 
one -dimensional sequential development which has been adopted will provide a direct 
access to the fascinating world of modern space dynamics. Starting from the most 
basic inverse -square law, the law of gravity for a sphere is developed. Then the sub- 
ject of motion of point masses under the influence of such a sphere is taken up. The 
first chapter thus serves as a necessary background and covers what is essential to 
the understanding of the rest of the book in the way of classical orbital theory. The 
development format used here is thought to be important in that notation is developed 
which will be used with a fair degree of consistency throughout the rest of the book. 
In addition, the groundwork for the methodology which will be adopted in succeeding 
chapters is laid. This methodology consists of a development from the Lagrangian to 
the resultant equations of motion to the solutions by way of whatever approximations 
are necessary. It is felt that this is the most direct and error -free way to go. In the 
event that the Lagrangian formalism may not be familiar to some readers, a suitable 
starting point is possible at the equations of motion which generally follow the 
Lagrangian within a few steps. There has been no attempt to place the book in a math- 
ematical straitjacket. Vectors and vector operations are used whenever they seem 



to be most appropriate, and thus clearest to the average student of mechanics. 
Emphasis has been placed on an intuitive feel for the concepts of energy and angular 
momentum which, it is hoped, most readers have acquired by the time they read this 
book. 

The second chapter is devoted to reentry theory and the relationship between 
reentry theory and the theory of the orbit. An attempt is made to move the reader 
away from classical passive observational theory to the developing modern ideas of 
active spaceships and controlled maneuvers. Emphasis is placed on the character of 
the gradual transition from pure space flight to pure atmospheric flight. Again, the 
reader should be constantly aware of the overriding importance of the ideas of angu- 
lar momentum and energy. 

Chapters three and four are concerned with relative motion between two bodies in 
orbit about the same planet, usually the earth or moon. In chapter three moving rec- 
tangular coordinate systems are taken up for the first time. This development leads 
naturally to the equations of relative motion between two simple point masses and, 
hence, to the Clohessy -Wiltshire rendezvous equations. The rest of this chapter is 
devoted to further attempts to improve these equations by the consideration of higher - 
order terms. In this chapter, perturbation solutions are discussed with emphasis on 
the physical content rather than on pure mathematical theory. 

The main concern of chapter four is also with relative -motion equations. The 
coordinate systems under consideration, however, are nonrectangular. It will be 
seen that most of the equations which are expressed in chapter three can be improved 
considerably in accuracy, or elegance, or both, when the coordinate system is chosen 
in a way that is more natural to the physical systems under consideration. The first 
of these systems is spherical and represents very little improvement except for the 
utility of the approach for certain pieces of hardware. The second system is in shell 
coordinates, which is essentially a set of cylindrical coordinates with the origin 
located at some point radially removed from the usual origin. It will be found that 
this system produces the best long-time or large distance -of -separation accuracy 
because terms of the series expansions can be retained to higher -order accuracy 
more readily and handled more easily. Again Clohessy -Wiltshire type rendezvous 
equations are developed and, in addition, the reasons for physical breakdown of 
relative -motion equations are carefully explored. 

Chapter five reverts to rectangular coordinates and treats the mechanics of 
simple rigid bodies under the influence of a gravity -gradient field, that is, the differ- 
ence in gravity between two points in close proximity. By the appropriate use of con- 
straints, this technique affords an introduction to rigid-body dynamics for satellites. 
A simple dumbbell is analyzed, and then the theory is extended to include the Interac- 
tions of two dumbbells. This approach allows a logical development to simple rigid- 
body dynamics including gyroscopic moments. Thus the book essentially covers the 
ground between point-mass mechanics for a body in orbit and rigid-body dynamics 
for space vehicles. At this stage, the possible directions of research are open-ended, 
and a development in terms of moments of inertia, a field which forms a separate 

vi 



subject for study in itself, becomes natural. With this subject the book is ended. It 
is emphasized that this is only a beginning and that almost every topic covered in 
this book affords a jumping -off place for more detailed subsequent research. It is 
hoped that this book will afford a background for continued improvement in the theory 
of space mechanics. 
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The Apollo 9 Command/Service Module as photographed in orbit above the earth. 



CHAPTER 1 



TWO -BODY MOTION 



Inertial or Newtonian Reference Frame 
Since to a large extent, the subject matter of this book is devoted to the study of 
orbital me'chanics in a variety of reference frames, the idea of an inertial reference 
frame will * a central concept. An inertial reference frame was once thought to be 
one that is fixed relative to absolute space. With the advent of relativistic mechanics 
the idea of absolute space lost its meaning, and the idea of an inertial frame of refer- 
ence was shown to be a property peculiar to an individual observer. Attempts were 
made to define a reference frame which is fixed relative to every body in the universe 
as an inertial reference frame. This idea is also not satisfactory, however, since 
every body in the universe seems to be in motion relative to every other body in the 
universe both in space and time. In addition, the interrelationship of mass and space 
is not satisfactorily understood at the present time. In the face of this vague founda- 
tion, an inertial reference frame will be defined as one in which Newton's laws of 
motion hold. It may be that in reality no such frame exists. There are situations, 
however, which approximate such a system. One is a reference frame fixed relative 
to the center of mass of our galaxy. For a large number of problems, including all of 
those which will be of interest in this book, a coordinate system fixed relative to the 
solar system or fixed relative to the center of the earth will be taken as a sufficient 
approximation to an inertial system. 

Assume then that an inertial frame x,y,z exists and suppose it is desired to 
compute the potential between a unit point particle and a finite rigid body as shown in 
the following figure: 




Figure 1-1.- Potential of a finite rigid tody. 
The potential between the point particles dM and m is defined as 
/ x n f dM(x',y',z') 

where G is some as yet unknown constant. 



(1-1) 
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Figure 1-2.- Infinitesimal mass unit and coordinates for computing the 
potential of a test mass m due to the presence of a sphere. 

Suppose now that the body M is a sphere of radius a and uniform density as 
shown in figure 1-2. Since the sphere is symmetric with respect to both spherical 
coordinates 6 and <p, the potential at any point a distance r from the sphere is 
given by equation (1-1) with 



dM = ap 2 sin dcp d(9 dp 
r2 = p2 + r 2 _ 2pr cos i 



(1-2) 
(1-3) 



where p is the radius over which the integration is to be performed and a is the 
density of the material of the sphere. Then 



u = G 



M 



dM 
R 



>a nir nZir 



G 



qp^sin 9 dcp d6 dp 



-'0 *^0 ^0 (P 2 + r - 2pr cos 9) 



1/2 



(1-4) 



= 2ttoG 



a n 



p 2 sin 6 d6 dp 



Jq (p 2 + r 2 - 2pr cos 6) 



1/2 



In order to facilitate integration at this point, let 

t = cos 6 

dt = -sin 9 d0 



(1-5) 



(1-6) 
(1-7) 



Then, 



u-^oopr 1 , ^«p 1/2 

J J t (p 2 + r 2 - 2prt) 



= ~2noG 



a 

P 



(p 2 + r 2 - 2rp) - (p 2 + r 2 + 2rp) 







dp 



2 i^J a p[(r-p)-(r + p)]dp 



It can therefore be seen that the potential of a sphere of uniform density is the same 
as that of a point particle of mass M at the center of the sphere. Thus, spherical 
planets may be replaced by mass points from the standpoint of dynamics. This pro- 
cedure will be followed for the most part in subsequent development. 

Newton's Law of Gravity 

It may be noted that the resultant force per unit mass between M and m is 
given by 



E = Vu=-GM(L 
m r 2 \r 



acting along a line joining m to M, or 



F =.9E^L(£j (1.8) 



r 

where G will now be defined as the universal gravitational constant. This is, of 
course, the law of universal gravitation first formulated by Newton: 

"Every particle of finite mass in the universe attracts every other particle with 
a force acting along the line joining them with a magnitude directly proportional to the 
product of their masses and inversely proportional to the square of the distances of 
separation." 

This law could be taken as a starting point for a theory of orbital mechanics. 
It was decided to take the potential approach, however, out of deference to field-theory 
concepts, which seem to be in the ascendency throughout other branches of mechanics 
at the present time. 

The classical two-body problem is developed next. The starting point for this 
development will be by means of the Lagrangian of a mass m in a potential field. 



Lagrangian of a Particle in a Gravity Field and 
Two -Body Equations of Motion 
The general expression for the Lagrangian in a noninertial frame of reference 
without translation is given by 

L = i mV 2 + mV • fi X F + I m( fi X F) - -t- (1-9) 

2 2 j r | 

wher e -^- is defined as the potential energy, V is the inertial velocity of m with 

I r*l — 

respect to the center of the potential field (i, and fi is the rotation vector of the 

frame of reference. By assuming, for the present, that Q, = 0, and hence the coordi- 
nate frame is nonrotating, the above equation reduces to the simpler form 

L l m y 2 _JL. (1-10) 

2 I r I 

In a cylindrical coordinate system, with out-of -plane coordinate z, 

V = Tr + ~fr6 + kz (1-11) 

1/2 
where the out-of -plane magnitude of r is given by |r|=(r2 + z 2 ) . The gravi- 
tational potential is given by the sum of the potential fields of both m and M. 
Hence, in this special case, 

H _ (MG. + ™G\ (1-12) 

or 

H _ G(m + M) (l-12a) 

i ? r ( r 2 +z 2) i/2 

where M is the mass of the central body and m is the mass of a satellite. Equa- 
tion (1-10) becomes 

L = I m(r 2 + vH 2 + z 2 ) + G * m + ^ (1-13) 

(r 2 + z 2 ) 

The equations of motion which result from this Lagrangian are 

„2 2 , G(m + M)r _ n (1-14) 

mr - mrfl + — * 1/2 

(r 2 + z2) 
mr 2 + 2mrr0 = (1-15) 

(r2 + z 2 ) 3/2 

Since G, m, and M are all constant, the gravitational potential at the surface of a 

planet can be defined by 

mg - r - G ^ M + m ) (i- 17 ) 

e e e r„ 



where r e is the radius of the planet and g e is a derived quantity known as the local 
gravitational acceleration. Thus, mg g r e is the potential energy at a distance r e 
from the origin. If g e is defined in this manner, 

G(M + m) = mg e r e 2 (1-18) 

If the right-hand side of equation (1-18) is now substituted for the gravitational term in 
equations (1-14), (1-15), and (1-16), the equations of motion become 

9 m g D r p 2r 
fnr-mrr + e Q = (1-19) 



(r 2 + z 2 ) 



3/2 



o.. . • 

mv^e + 2mrr0 = (1-20) 

mg r 2 z 
mz+ S ^2 = (1-21) 

(r 2 + z 2 ) 

These equations completely define two-body orbital motion in an inertial frame of ref- 
erence. It should be noted that equation (1-20) is actually a torque -balance equation 
while equations (1-19) and (1-21) are force-balance equations. Equation (1-19) only 
becomes a force equation when divided by r. The peculiar nature of equation (1-20) 
is readily seen to be a consequence of the choice of a cylindrical coordinate system. 

Integration of equation (1-20) gives the angular momentum associated with the 
plane containing r and 9: 

m^)=0 (1 _22) 

Let 

r 2 = h (1-23) 

where h is called the angular momentum per unit mass. From a strict vector 
viewpoint h is only one of the components of angular momentum, that is, that com- 
ponent formed in the direction of TxJ = k in the notation of equation (1-11). Since 
this component is conserved irrespective of the other two, however, it is especially 
valuable in orbital mechanics. 

Equations (1-19) and (1-20) can be considerably simplified by assuming that z 
and z are initially zero. If this assumption is made, equation (1-21) remains zero 
for all time and the motion remains planar. Equation (1-19) becomes 

•2 m S P r e 2 
mr - mrr + — x. — =o (1-24) 

T l 

Equation (1-20) is not affected by neglecting z and z; so equation (1-23) still 
applies. Under these conditions there is no angular momentum in the z equation 
and, since r is radial and so cannot express angular momentum, all of the angular 
momentum is contained in h. Since it is always possible to select a coordinate 



system where the out-of -plane motion is zero initially, all simple orbits can be treated 
as planar so long as there is no acceleration acting in the out-of -plane direction. 

Substituting the right-hand side of equation (1-23) into equation (1-24) uncouples 
the latter, thus giving 

h 2 mg r e 2 

mr - m -r + - — >- — =0 u-zoj 

This equation can be multiplied through by r and integrated to yield 

1 o 1 u2 mg r p * 

I mi . 2 + I m h---^-l-= Constant (1-26) 

2 2 r 2 r 

Let the integration constant be designated by mH. The left side is the law of conser- 
vation of energy, so H is the energy per unit mass and mH is the total energy. 
Thus 

i mr 2 + I m 4 . mg e r e 2 „ mH (1 . 27) 

2 2 r 2 r 

The orbital speed at any given radius may be found from equation (1-27). If V 
is speed, 

ImV 2 =im[r 2 +(re) 2 ] (1-28) 

Then using the left-hand side of equation (1-28) in equation (1-27) and dropping the 

common m 

2 

1 V 2 = H + ^- (l- 2 9) 

2 r 

for all orbits. 

By setting r = the existence of extremal distances can be investigated. 

Equation (1-27) is solved for r to give 



-l ± Jl + . 2h2H 



r, 



2H/g p r f 2 



e e (1-30) 



3 e e 



It can be seen, since g e and r e are always positive and h is even-powered, that 
the sign of H governs the nature of the orbit in a distinctive manner. If H is 
zero or positive, there can be only one positive nonzero extremal. If H is negative, 
there are two real values. In this case the positive sign of the radical gives an 
expression for minimum distance from the origin r p , the perigee radius, while the 
negative sign gives the maximum distance r a , the apogee radius. 



Equation of a Conic Section 
Equation (1-27) can be integrated with a change of independent variable from t 



to 6 to obtain the classic equation of a conic section. By equation (1-23) 
A h 



So 



r = 



h_dr 
r 2 d6 



(1-23) 



(1-31) 



Equation (1-27) becomes, after some rearrangement, 



drr = r 2 



2 S^ e 



Taking the square root, the variables separate to give 



(1-32) 



de = 



dr 



'^\r2 + 



2g_r 



e e 



r-1 



1/2 



(1-33) 



This equation can be integrated directly (ref. 1) to give 



= cos 



r 5 e e 







; e 2r e 4 + 2h2 H ) 



1/2 



where e o is an integration constant. If this equation is solved for 



(1-34) 



r, one obtains 



r = 



h2 /Ve 2 



1 + 11 + 



2h^H 

o- 2 r 4 i 
°e e / 



.1/2 



cos 



(fl - e ) 



Setting 



p = 



g e V 



which is called the semilatus rectum and, if e is defined by 



(1-35) 



(1-36) 



6=1 + 



9 \l/2 
2h^H \ 



g e Vy 



(1-37) 



equation (1-35) becomes 
P 



r = 



1 + e cos(6 - O ) 



(1-38) 



Examples of typical solutions of equation (1-38) with p set equal to one are shown 
in figure 1-3. The same information in more familiar rectangular coordinates is 
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Figure 1-3.- Typical conic sections in nondimensional form. Radial distance 
plotted against angle in degrees. 



shown in figure 1-4. Equation (1-38) is a standard form of a conic section with the 
origin at one focus and eccentricity e. The three types of orbit suggested by equa- 
tion (1-30) can now be classified by the value of the eccentricity e: 

If i e < 1, the orbit is an ellipse 

If e = 1 , the orbit is a parabola 

If e > 1 , the orbit is a hyperbola 
In the special case where e = equation (1-38) reduces to the equation of a.circle. 

These cases will now be examined in detail with special consideration given to 
the elliptic case since this is the most important case from an astrodynamic point of 
view. 
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-3.0 -2. 



p= I all cases 




-25 



Figure 1-4.- Conic sections in rectangular coordinates. 
Elliptic orbit .- In the first case above 

vl/2 



< e = 1 + 



2h 2 H 



< 1 



or 



2h 2 H 
>e e 



< 



(1-39) 



(1-40) 



Since h, g e> and r e are all even-powered, H must be negative; that is, the total 
energy per unit mass (with respect to a zero potential infinite point) is negative for an 
elliptic orbit, and the spacecraft is trapped in a potential well. 

11 



The next step is to calculate the period of an elliptic orbit. This period is fre- 
quently obtained by introducing coordinates associated with the geometric center of 
the ellipse and making an appropriate transformation. The period is then found in 
terms of the semimajor axis distance. A more direct, though algebraically more 
involved, method is to calculate the area of the ellipse divided by the rate at which 
area is swept out by the r vector. In this instance the orbital period P is 



1 



2tt 



de 



dt = \ =r = 
J ^0 r*6 



271 r 2 d0 



2tt Jl 



?ld£ 



(1-41) 



The numerator of the above expression is found by multiplying equation (1-33) by r^ 



r 2 de = • 



r dr 



2k r„ 



^M^?- 1 



1/2 



Hence, taking advantage of symmetry over r in the limits 



(1-42) 



2tt Jl 



= J "h 



d6>_2 
h 



r dr 




w^-\ 


1/2 



(1-43) 



where r a is apogee and r p is perigee. 

It can be seen that by applying the limits in the above equations a tacit assump- 
tion is made that the orbit is closed and hence is an ellipse. Integrating from tables 



p - h. 
H 



r 2 + 



2 g fi r , 



e e 



1/2 



g e r e . -1 

+ _£___ sin 



HV-2H 



2Hr + g e r e ' 



(2Hh 2 + g e 2 r e 4 ) 



1/2 



(1-44) 



H 



fc*m*-> 



1/2 



2H\ 2 , 
2 P 



'2g r 



e e 



1/2 



S e r e . -1 

+ : sin 



H\P2H 



2Hr a + gr e ' 



(2Hh 2 + g e 2 r e 4 ) 



g p r e -1 

e e sin L 



nfM 



2Hr p + geV 



(2Hh 2 + g e 2 r e 4 ) 



1/2 



(1-45) 



Then, making use of the values of r a and r p in equation (1-30) the first two terms 
in equation (1-45) vanish and the second two terms reduce to 



12 



p = 



g e r e 2 -1, g e r p 2 1 g D T Jn 

_^-£-. sin x (l) + -£-£- sin X (-l) = - e e 



HV^2H 



HV-2H 



HV-2H 



(1-46) 



A negative period, of course, corresponds to a reversal of direction or of time; hence 
only the absolute value is of physical interest. This equation brings out the fact that 
the period of an orbit depends only upon the total energy of the orbiting body and the 
strength of the gravity field of the attracting planet. 

The semimajor axis of an ellipse a is given by one -half the sum of the perigee 
radius and the apogee radius. Again making use of equation (1-30) 



r a + r p 



, 1 f -ge r e 2 - \/geV + 2Hh2 | -g e r e 2 + feV + 2Hh2 ) 
2 \ 2H + 2H / 



2 2 \ 2H 2H 

The period in terms of the semimajor axis becomes 



£ r * 
5 e e 

2H 



(1-47) 



P = 



2ir 



-g_r_2 



e e 



2ira 



= 2ir\ 



\f-2H\ 2H / ^h "Vg D r D 2 

Ice 



(1-48) 



It can be seen that the period becomes a function of either (a,H) or (a,g e r e 2 ). About 
a given planet, any two orbits having the same energy and semimajor axis are syn- 
chronous. This will be found to be important in the study of rendezvous problems. 
From equation (1-36) it can be seen by rewriting the equation that 



P = 



g r ^ 

& e e 



2H 




g r < 
& e e 

2H 



(l - e 2 ) 



(l-36a) 



And so, from equation (1-47) 

p = a(l - e 2 ) 
Thus equation (1-38) can be written for an elliptic orbit in the useful form 
a(l - e 2 ) 



r = • 



(l-38a) 



1 + e cos( 6 - 6 Q \ 
Circular orbit .- A special case of an elliptic orbit is a circular orbit. In this 



kind of orbit e = and equation (1-35) reduces to 



cir 



3 e e 



■ = Constant 



H = 



5 e e 
2h 2 



ge r e 
2r cir 



(1-49) 



The energy of a circular orbit can be found from equation (1-37). Since e is zero 

2 



(1-50) 
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The speed is found very simply from equation (1-29) and is 

2 



I v 2 =- ger& 



2 cir 2 r cir 



or 



V • 

cir 



g e r e (1-51) 



r . 
cir 



The angular momentum can also be established from equation (1-50) and is given by 



u - r \fek_-r ■ V (I- 52 ) 

n cir - r cirU r cir " cir cir 

Hence it can be seen that 9 is constant. Denoting this constant by w gives 



h . =r 2. w = r . \!£e! (1-53) 

"cir x cir cirU r cir 



Hence 



v 6 . 
cir 



g e r e (1-54) 



All of the other important relations for circular orbits can be derived in a similar 
manner from the appropriate elliptic equations. 

Parabolic trajectories .- For parabolic trajectories e = 1; therefore, from 
equation (1-37) 

2h 2 H _ (1-55) 

or 2 r 4 

s e e 
and from equation (1-35) 



h 2 /g e r ( 



2 
e 



(1-56) 



l + cos(e - e ) 

Since h 2 and g 2 r 4 are finite, and nonzero, the total energy, from equa- 
tion (1-55), is zero; thus 

H = o d- 57 ) 

The zero total-energy state could be considered the main distinguishing feature of a 
parabolic orbit. The speed is found again from equation (1-29) and is e 

1 v 2 ■■■. geFe 

2 parabolic r 
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or 



V, 



2g„r 



ic-f-V 1 " (!- 58 ) 



parabolic 

A simple and useful relation can be found between circular and parabolic speeds by 
comparing equations (1-51) and (1-58) for trajectories at the same altitude. If 

r = r cir' 

* Vrabolic = ^ (15g) 

^circular 

The minimum distance from the origin for a parabolic trajectory occurs when r = 0. 
From equation (1-27) by setting the time derivative of r equal to zero 

h 2 

'min^r^ (1 " 60) 

2g e r e 2 

Thus, if one wishes to inject a vehicle from circular orbit to a parabolic trajectory by 
thrusting tangent to the orbit, the required angular momentum is 



\ 2g e r e 2 
h parabolic ~ r min\rJ : (1-61) 

I mm 

The ratio of angular momenta can also be found. From equations (1-53) and (1-61), 

Vrabolic g ^ (1 _ 62) 

circular 

The difference in angular momentum between h parabolic and h circular is given by 
slight rearrangement, making use of equation (1-52) 



& - btrff' 



parabolic " Circular _ ^ 2 " 1 J r cirU~F~ ( 1_63 ) 



A parabolic trajectory has the property that a body on such a trajectory has just 
barely enough speed to escape from the gravitational field of the attracting planet. 
Thus these trajectories are often called minimum -escape trajectories. Such trajec- 
tories are not encountered very often in actual practice. However, bodies entering 
the solar system with essentially zero velocity relative to the sun are found to approx- 
imate this situation fairly well. Included among these are a large class of comets. 
Since these bodies return to the same distance from which they start, comets in per- 
fectly parabolic orbits are observed only once. Due to minor energy losses, however, 
some of these bodies are trapped by the solar system and thereafter assume highly 
elongated elliptic orbits. Under these conditions, they periodically return to the 
vicinity of the sun. 

Many meteorites are also presumably on parabolic trajectories with respect to 
the earth. They travel around the sun in essentially the same orbit as the earth until 
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they encounter the earth's gravity field. In this case, few of the observed ones sur- 
vive an encounter with the earth. 

Hyperbolic trajectories .- Hyperbolic trajectories are important in interplane- 
tary work because bodies in such orbits travel with speeds in excess of minimum - 
escape speed. Up until the present time no such trajectories have been attempted 
relative to the solar system owing to the high energies involved with respect to the 
energy level of the earth. Such trajectories are important in interplanetary work, 
however, since they allow the shortest travel times between planets. Of course, any 
interplanetary probe leaves the vicinity of the earth along a hyperbolic trajectory 
with respect to the earth itself. 

For a hyperbolic trajectory e > 1; therefore, 

/ 2 \ 1/2 

/ 1+ _2hfH_\ >x (1 _ 64 ) 

Since, as before, h, g e , and r e are all even-powered, H must be positive; that 
is, the total energy per unit mass (with respect to a zero potential infinite point) is 
positive, and a spacecraft has the capability of leaving the local gravitational potential 
well. A body having such energy with respect to a planet's gravitational field escapes 
into the solar system; a body having such potential energy with respect to the sun's 
gravitational field escapes from the solar system into interstellar space. 

Equation (1-35) applies as in the case of an ellipse although the results are, of 
course, quite different, as zero denominators are possible. Thus 



'g r 
r = 



h 2 /| 



2 
e ~ e (1-35) 



1+ rfv) cos(e "'' o) 

The concept of period is meaningless, as such orbits are not cyclic. 

Mean and Eccentric Anomaly 
The mean motion of a satellite in elliptic orbit is defined as 



n _ 277 = 2HV3JH (1 . 65) 

P _cr r " 

s e e 
or, expressed in terms of the semimajor axis, from (1-47) 



n = 




(1-66) 



This is the rate at which an orbiting spacecraft would have to move if its angular 
velocity were constant in order to complete one orbit in time P. The mean angle 
traveled from some arbitrary starting time t would then be 
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2H\PZH'(t - t 



M = n[t - t, 



i l ~ l o) 



f2H(t - t Q ) 



-E r ^ 
& e e. 



(1-67) 



The eccentric anomaly E, described in figure 1-5, is defined by the transcen- 
dental equation 



E - e sin E = M = 2H ^' 2 ^ (t - t ) 
& e e 



(1-68) 



or, in terms of the semimajor axis, 



E - e sin 



E = M = v e e V '- = n(t - t„\ 

a 3/2 I °) 



(1-69) 



This equation is called Kepler's equation. 



Projected circle 



Actual orbit e = 0.5 




Vernal equinox, T 



Figure 1-5.- Eccentric anomaly description. 
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Equation of the Center 

Equation (1-69) can be used to derive an important series expansion of an 
ellipse called the equation of the center which will be needed in chapter 3. The par- 
tial derivative of equation (1-69) with respect to the eccentricity is given by 



8E . r, „ 8E 8M 

— - sin E - e cos E — = — — 
de 9e 8e 

From the second equality in equation (1-69), it can be seen that 
Bift - t > 



(1-70) 



8M _ g_ 
de 8e 



= 



(1-71) 



so that 



or 



8E 



(1 - e cos E)— - sinE = 



8E 
8e 



8e 



e sin E 
1-6 cos E 



Also, from equation (1-69) the exact differential is given by 

dM = (1 - e cos E)dE 
Equations (1-72) and (1-73) can be combined to give 

e il dM = e sin E dE 

de 

This equation can be integrated over M to give 
>M 



e f ££. dM = -e cos E + C 
Jn 



8E 
)q de 



(1-72) 



(1-73) 



(1-74) 



(1-75) 



where C is an integration constant. In order to perform this integration, it is nec- 
essary to expand E in a series of terms in M. To facilitate the integration to suf- 
ficient accuracy, solve equation (1-69) in a nested expression: 



E = M + e sin E 



but 



sin E = sin(M + e sin E) 

= sin[~M + e sin(M + e sin E)1 
= sinj~M + e sin(M + e sin((. . .)))] 
Then apply the trigonometric identity for the sine of the sum of two angles 
sin(A + B) = sin A cos B + cos A sin B 



(1-76) 
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to obtain 

E = M + e sin MJ~cos(e sin(M + e sin(. . .))f) 

+ e cos M[sin(e sin(M + e sin(. . .)))] (1-77) 

But for small eccentricity this equation becomes, approximately (ref. 2) 

E ~ M + e sin M + e 2 cos M sin M + Higher-order terms 
so 

* e 2 

E ~ M + e sin M + — sin 2M + Higher-order terms (1-78) 

The partial derivative of equation (1-78) with respect to e can then be applied in 
equation (1-75), resulting in 

-e cos E = C + e \ (sin M + e sin 2M + . . .)dM (1-79) 

J 

where C is still an undetermined constant. From figure 1-5 

ae + r cos(6> - S ) = a cos E (1-80) 

and from equation (l-38a) 

a(l - e 2 ) 
r = - V L's (l-38a) 

1 + e cos (0 - 6 ) ' 

If (9 - 6 ) is eliminated between these two equations, the equation of an ellipse in 
polar form is given by 

r = a(l - e cos E) (1-81) 

Hence 

-6 cos E = - - 1 
a 

Substituting this value for the term on the left in equation (1-79) 

J 

or, integrating, 



i = 1 + C + e \ (sin M + e sin 2M + . . .)dM (1-82) 



\= 1 +C - e cos M -ie 2 cos 2M (1-83) 

If it is assumed that the series is uniformly convergent for small eccentricity, this 
equation can again be integrated between the limits M = and M = 2tt: 

f 2tt r 2ir 2 p27T 

\ - dM = 2tt(1 + C) - e \ cos M dM - 1- \ cos(2M)dM (1-84) 

J J 2 J 

The two right-hand trigonometric terms vanish and 

f 27r r 

| dM = 2tt(1 + C) (1 - 85) 

J 
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The left-hand side is evaluated by reverting back to eccentric anomaly form, noting 
that the limit remains unchanged, 

f V - € cos E) 2 dE = 2ir(l + |-j / (1-86) 

Thus 

2tt(i +~) = 27r(l+C) 

r- e2 (1-87) 

I=l-ecosM-i e 2 (cos 2M - 1) + Higher-order terms (1-88) 

a 2 

The computation of the appropriate series for 6 follows from this equation and 
equations (l-38a) and (1-69). By eliminating r between (1-88) and (l-38a) 

l-ecosE = - l -^4—ir\ (1 - 89) 

1 + € cos(0 - B ) 

This equation can be solved for cos(e - 6 ) by rearrangement 



and 



Then 



de = 



(l - e2) 1/2 



cos(e _e o) = J^JLzJL. (1 . 90 ) 

v ; 1 - e cos E 

and can also be converted to sine form, which is 

/ o\l/2 

(1 - e^J sin E ,* Q1 x 

sin 6 - e ) = *-- '- — — U-91) 

v °> 1 - e cos E 

Equation (1-90) can then be differentiated with respect to E 

SAa (l-e 2 )sinEdE (ig2) 

sinfe - e )de = ■* - — (i~ 92 ) 

(1 - e cos E) z 
Elimination of sin(6 - O ) between equations (1-91) and (1-92) results in 

1 /2 
dc (l - e 2 ) (1 - e cos E)dE (^93) 

(1 - e cos E) 2 
which is, from Kepler's equation, equation (1-73), 



■dM (1-94) 



(1 - e cos E) 2 
This equation becomes, as a consequence of Kepler's equation, r 

de=\tT7W 2 dM (1-95) 

IdM/ 
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The quantity dE/dM is found at once from the expression for E in terms of M, 
equation (1-78), by taking the derivative of E with respect to M, 

dd = (l - e 2 ) (l + e cos M + e 2 cos 2M) dM (1-96) 

/ 9\!/2 
Then, after expanding ^1 - e*) in the power series 

(l - e 2 ) _ i _ L- + Higher-order terms (1-97) 



2 



de = 1 1 - — 



-U-4?: 



) - 2 - 

(l + e cos M + e 2 cos 2MJ + Higher-order terms 

J 1 + 2e cos M + e 2 (2 cos 2M + i- cos 2M + -) + Higher-order terms 



dM 



1 + 2e cos M + 1 e 2 cos 2MjdM (1-98) 

Integrating with respect to M, 

5 9 

9 = M + 2e sin M + - e^sin 2M + Higher-order terms (1-99) 

Equations (1-84) and (1-99) describe the motion of a body in an elliptic orbit in terms 
of the mean anomaly. These equations diverge for eccentricities greater than 0.6627 
and converge rapidly for typical orbital eccentricities and hence are handy for manual 
calculation. They will be used for analytical expressions for noncircular orbits in 
chapter 3. Equation (1-99) is referred to as the equation-of -the -center in some books 
on orbital mechanics. 

The Orbit in Space 

Discussion has thus far been confined to that of a plane two-dimensional orbit 
such as, for example, an ellipse. It can be seen that the ellipse is completely speci- 
fied in its plane by four independent, fixed parameters; for instance, from equa- 
tion (1-35) the conic section could be completely specified by naming the angular 
momentum h, the energy H, the gravity parameter g e r e 2 , and the argument of 
perigee O . 

The orbit could just as well have been defined by the following parameters: the 
length of the semimaj or axis a, the eccentricity e , the argument of perigee 9 , and 
the time of perigee passage. 

In each of these instances it is found that there are three physical parameters of 
an ellipse plus one quantity which specifies the orientation of the ellipse in the refer- 
ence plane. In addition to these, two quantities are required to specify the orientation 
of the plane of the orbit with respect to three-dimensional space. These quantities 
are usually specified as the inclination, which is the angle the plane of the orbit makes 
with respect to the fixed reference plane; and the longitude of the ascending node, 
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which marks the orientation of the line of intersection of the orbit with respect to the 
reference plane. There are a variety of such systems in existence, each having its 
advantages for certain problems. Following the main line of reasoning used here, 
however, the system which falls naturally out of classical mechanics will be used; 
namely, energy, angular momentum, gravity, and whatever geometric orientation 
coordinates seem most appropriate for a given problem. Orientation can be assumed 
to be specified by inclination and the longitude of the ascending node. Under this 
organization form, a coordinate system such as is shown in figure 1-6 is appropriate, 
where 9 is the argument of perigee, i is the inclination of orbit plane to plane of 
the ecliptic, and fl is the longitude of the ascending node. 

y 



Line of nodes 




Plane of motion 
of the satellite 



Equatorial 
reference plane 



Ascending 
node 

Figure 1-6.- Orbital elements in space. 
Assume a rectangular coordinate system with the x-, z-axes located arbitrarily 
in the equatorial plane and the y-axis perpendicular to it. Then the position of the 
satellite at any given time is specified by 

x = rfcos(0 - 6 )cos fi + sin(0 - # )cos i sin fij 

y = r[sin(e - o )sini] \ (1-100) 

z = r|sin(0 - o )cos i cos O - cos(6 - 6 ) sin ^j 

where 9 is the true anomaly and r is the radial distance to the satellite. 

Since 6 , i, and Q, are all invariant for a planar orbit, it remains only to 
specify the relationship between r and 9 to define completely the position of the 
satellite. This relationship can be specified by the angular momentum or the 
Lagrangian. 
22 
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CHAPTER 2 



REENTRY THEORY 



The objective of this chapter is to express the equations of orbital motion in a 
coordinate system which is natural for a reentry vehicle and to make an attempt at a 
solution. In so doing it will be necessary to introduce the concepts of aerodynamic 
lift and drag. First, however, it will be found convenient to express the orbital 
motion in terms of the velocity vector of an orbital vehicle. Efforts shall be confined 
to motion in the orbital plane, keeping in mind that in actual practice out -of -plane 
motion can be obtained either by thrusting out of plane or, in the atmosphere, by roll- 
ing the entry vehicle so as to produce an out-of -plane lift component. In the orbital 
plane the velocity vector is specified completely by the magnitude V and an orienta- 
tion angle with respect to the local horizontal called the flight-path angle y. 

It is easy to see why this formalism has been useful over the years. While in 
this form it is true that radial velocity and angular momentum are not uncoupled, the 
drag force is directed opposite to that of the velocity vector. The magnitude of this 
force is found to be approximately proportional to the projected surface and the 
dynamic pressure acting on this surface. Thus, as can be seen in the following 
sketch (fig. 2-1), a simple resolution of forces in the direction of motion is possible. 



Lift 




Axis of reentry- 
vehicle 

Figure 2-1.- Lift and drag forces on a reentry body. 

The Lagrangian for a body in orbit about the earth in an earth-centered, iner- 
tial, polar coordinate system (z neglected) is 

L= l m (r 2 + A 2 ) + ^1 (2-1) 

The equations of motion which result are # 

2 mg g r 2 
mr - mrfl + — = (2-2) 



mr^'S + 2mrr0 = 



(2-3) 



24 



Note that equation (2-2) is a force equation while equation (2-3) is a torque equation 
and only becomes a force equation when divided through by r. These equations are, 
of course, the same as equations (1-14) and (1-15) of chapter 1. 

Speed and Flight-Path Angle 
Now make use of the relations which define speed and flight-path angle 

M V cos y = r0 (2-4) 

V sin y = r 



(2-5) 



It follows from equations (2-4) and (2-5), by taking the time derivative, 



■A • • • • • V 2 

r 9 = V cos y - Vy sin y - re = V cos y - Vy sin y - -L- sin y cos y (2-6) 



Hence equation (2-3) becomes 



mV cos y - mVy sin y - m ^—- sin y cos y + 2m ¥— sin y cos y = 



or 

r2 



mV cos y - mVy sin y + m ■— sin y cos y = (2-7) 

Also, making use of equations (2-4) and (2-5), equation (2-2) becomes 

o p" -)-» & 

mV sin y + mVy cos y - m ^— cos y + m e e = (2-8) 



2 

Now it is necessary to combine equations (2-7) and (2-8) in such a way as to separate 
V and y: 

(1) Multiply equation (2-7) by cos y and equation (2-8) by sin y 

mV cos y - mVy sin y cos y + m -^- sin y cos 2 y = (2-9) 



mV sin y + mVy sin y cos y - m %- sin y cos 2 y + m -—- sin y = (2-10) 



r^ 



,,„2 



Adding equations (2-9) and (2-10) and making use of sin y + cos y = 1, 

mV + m-^-siny=0 (2_n) 

(2) Multiply equation (2-7) by sin y and equation (2-8) by cos y 

mV sin y cos y - mVy sin y + m ^L S in 2 y cos y = (2-12) 






2 V^ 3 ^e e n in i«!i 

mV sin y cos y + mVy cos y - m -y- cos y + m — ^— cos y = U (i-ii) 



r 2 



2 2 

Subtracting equation (2-12) from equation (2-13) and making use of sin y + cos y = 1 

2 S r 

m Vy - m — cos y + m e e cos y = (2-14) 

or, since gr 2 = g e r e 2 for any r 

mVy - mg(pr - ljcos y = (2-15) 

Lift and Drag 

Equations (2-11) and (2-15) represent the motion in terms of velocity and flight- 
path angle of a body in orbital motion (refs. 3 and 4). It should be emphasized that 
this is no more than an algebraic transformation and as such contains exactly the 
same information as equations (2-2) and (2-3). All that is done is a conversion from 
orthogonal velocity components to a scalar magnitude and a direction angle. Although 
velocity flight -path-angle coordinates may be more intuitive than the formalism of 
equations (2-2) and (2-3), it should also be noted that there is a penalty which must be 
paid. The great beauty of the classical orbital equations, equations (2-2) and (2-3), is 
that the angular -momentum equation is uncoupled from the energy equation, thus per- 
mitting direct integration of the equations of motion. In velocity flight -path -angle 
form this situation is no longer true, and no exact differential for angular momentum 
exists. Equations (2-11) and (2-15) are each expressions of combinations of angular 
momentum and energy and have to be uncoupled by essentially reversing their deriva- 
tion. This fact can be seen very easily by examining the transformation equations, 
equations (2-4) and (2-5). It is clear that no expression equivalent to r § is pos- 
sible using just equation (2-11) or equation (2-15). Equations (2-14) and (2-15) 
express the reentry equations in a form which is applicable when there are no non- 
conservative forces present; that is, they are valid for space applications but are not 
valid in the atmosphere where lift and drag are present. It is thus necessary to add 
the nonconservative forces on the right-hand sides of these equations. The magnitude 
of the drag force is found to be approximately proportional to the product of the pro- 
jected surface and the dynamic pressure acting on this surface. Denoting the propor- 
tionality constant by Cr> then 

D = C D Sq ( 2 ~ 16 ) 

where the dynamic pressure is given by 

q = IpV 2 * ( 2 "17) 

2 

and p is the atmospheric density. The lift force is defined analogously by 

L = C L Sq (2- 18 ) 
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It should be noted that C L and C D are both only approximately constant but this 
approximation is found to be good experimentally. Hence, for orbital flight in an 
atmosphere, equations (2-11) and (2-15) along with equations (2-4) and (2-5) form the 

complete set 

mV + mg sin y = - |pC D SV 2 (2-19) 



mVy - mgfe - ijcos y = |pC L SV 2 



(2-20) 



with 



rd = V cos y 
r = V sin y 

2 



g = 



g r 
& e e 



(2-4) 
(2-5) 

(2-21) 



to specify the local gravitational attraction at the vehicle. 

The Barometric Equation 

In equations such as equations (2-19) and (2-20) it is necessary to have a rea- 
sonably valid expression for the density as a function of altitude. This expression is 
obtained by combining the hydrostatic equation with the perfect gas law. Consider a 
column of gas with the density known for every given altitude as shown in figure 2-2: 



2ZZ 



p 
22 



•P-(f )dh 



£ 



p 



77777777777/77777777777 



Surface of earth 



Figure 2-2.- Integration scheme to determine pressure 
as a function of altitude. 
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The hydrostatic equation is found by summing the forces acting on a small vol- 
ume of gas in equilibrium. From the figure 



P - (|j)dh - pg dh = p 



So 



lr- 



or 

dp = -pg dh (2-22) 

The hydrostatic pressure p for some altitude h is given by integrating the product 
of pressure and density from this altitude to infinity: 



p=\ pgdh (2-23) 

J h 

Equation (2-22) is to be combined with the perfect gas law. For any specified volume 
of gas which is assumed to be homogeneous and isotropic the perfect gas law gives 
the pressure as 

p = — p (Perfect gas law) (2-24) 

where 

k Boltzmann's constant, 1.380622 joules/K 

T absolute temperature 

M mass of air contained in unit volume 

Thus, if it is assumed for any sample volume that temperature and mass are constant, 
dp = ^dp (2-25) 

and equation (2-22) becomes 

|Idp=-pgdh (2-26) 

M 

Inserting the value of g from equation (2-21) and rearranging slightly 

2 



dp _ M g e r e 
P " kT r 2 



dh (2-27) 



Expressing r as 

r = r e +h € (2-28) 

and integrating from surface h = to some altitude h 



P d£- JL K r 2 C h dh 
p " kT 8 e e J 2 



(2-29) 
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Hence 

2 



In -£-= M £e r e 

P e " kT ( r e + h) 



(2-30) 
An equation is found by solving for p 



P = P e exp 



M g e r e 2 h 

kT (r e + h) r e 



(2-31) 



It is usually convenient at this point to make the approximation that r e * r + h 
since if is small compared to r e . Then a simpler form of the equation results 

/Mg 



where 



v^T/ h 
P = P e e (2-32) 

g r 2 
ga (r e +h)r e (2-33) 

Since the coefficient of h in equation (2-32) is very nearly constant, it is usual 
practice to treat it as a constant. Setting 

kf = ^ (2-34) 

the usual aerodynamic form of the equation becomes 

P = p e e '^ (2-35) 

where /3 is frequently expressed in reciprocal terms, 
J 
H 



= iT" (2-36) 

s 



The symbol H s in this form is called the scale height and is important for planetary 
atmosphere studies since it is the altitude at which atmospheric density decreases to 
1/e = 0.3679 of its value at the surface. For spacecraft entry, however, /3 is 
somewhat handier. For the earth the value of £ is 0.0001366/m and the scale height 
is approximately 7.3 km. 

The usefulness of the barometric equation is not so much in the numerical 
results it produces as in the qualitative insight that can be gained from it in the mech- 
anics of atmospheric entry. It can be inferred from the derivation that this equation 
does not yield highly accurate results, largely because the actual atmosphere is not 
isothermal as was assumed in the derivation. As a practical matter the equation is 
inaccurate by as much as a factor of two between 60 and 80 kilometers. This state of 
affairs is not as bad as it may seem, however, as it merely means that expected 
maneuvers will occur at times different from those expected and at different altitudes. 
Qualitatively an entry vehicle is not that sensitive to atmospheric-density variations 
for a large range of flight regimes. 
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One method of getting around this problem would be to divide the altitude of 
interest into layers and to determine values of /3 and p e which accurately describe 
the density variation for each small layer. The analytical solution to the equations of 
motion to be developed in the next section, which assumes an isothermal atmosphere, 
could then be used to solve each problem in a stepwise fashion. The more practical 
solution, however, when using a large computer is to use numerical density tables 
from the U.S. Standard Atmosphere, and this is what is usually done. 

Solution of the Atmospheric-Entry Equations 

Over the years a number of closed-form solutions of equations (2-4), (2-5), 
(2-19), and (2-20) have been published. All have the expected drawbacks of approxi- 
mate solutions to nonlinear differential equations in not being accurate enough for 
practical applications and in not being valid over a wide enough flight regime. These 
solutions have been of substantial value, however, in that they have provided a great 
deal of understanding of the nature of the problem. In addition, analytical solutions 
have been of some service in the development of guidance equations. Only two 
approximations are justified and they are 

(1) sin y ~ 

V 2 

(2) 1- - 1 « 

gr 
The first approximation (ref . 4) is an expression of the fact that most entries are 
made at a shallow flight -path angle in order to minimize heating and maximum accel- 
eration. While the total integrated heat at shallow flight -path angles is quite large, 
the maximum heat-input rate is small; hence such trajectories are the preferred ones 
from a heat -dissipation standpoint. The second approximation, though not necessarily 
limited in entry angle, is valid only wherever the velocity is equal to satellite velocity. 
A solution is presented here which is essentially the same as in reference 3, utilizing 
only the second assumption, primarily for the purpose of illustrating the general 
technique and for the sake of understanding. Solutions making use of assumption (1) 
above are abundant in the literature and none will be presented here. 

The atmospheric -entry equations can be written in the form 

mV = -D - W sin y 



mVy = L + wf^r - l) cos y 



€ 



) (2-37) 

h = V sin y 

r6 = V cos y 
where 6 is the down -range angle about the center of the earth and is included for 
the sake of completeness. A numerical integration of this set of equations for an 
unguided spacecraft having essentially the same physical characteristics as an 
Apollo spacecraft is shown in figure 2-3: 
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(b) Deceleration versus time. 
Figure 2-3.- Unguided Apollo reentry trajectory. 
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(d) Altitude versus time. 
Figure 2-3.- Continued. 
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(e) Total energy versus time. 
Figure 2-3-- Concluded. 

For simplification of the atmospheric -entry equations, the following substitu 
tions and assumptions are made: 

Assumption (1) 

The barometric equation for density applies 

-/3h 



where p and are constants, and 



(2-35) 



^=K e 



-0h v 2 



where K, 



C DPe S 

2m 



= Constant, and 



m"H e v 



CrOS 

where K<> = ■ ■ = Constant. 
6 2m 



Assumption (2) 



■—[—- ljcosy = e(V) ~0 
v 2\g r / 
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From the second assumption it is seen that whenever the velocity is equal to 
satellite velocity, that is, V = \fgr, or when y = 90°, the term e(V) is identically 
zero. Hence, these closed-iorm solutions to the reentry equations are developed by 
assuming that the difference between the centrifugal and gravitational accelerations is 
negligible. The term e(V) will be carried, however, until this assumption can be 
justified. 

Equations developed under these assumptions apply to a number of physically 
useful problems; for instance, these equations would apply to most entries which are 
steep entries in the uppermost regions of the atmosphere. They also would be appli- 
cable to skip trajectories for satellites and spacecraft and to grazing entries by such 
objects as meteors and tektites. 

For simplicity, the following change of the variable h to the variable y is 
made. Let 

y = e ~^ h (2-38) 

so that 

|=-/3h (2 -38a) 

With the foregoing assumptions and the change of variables given by equations (2-38) 
and (2 -38a), the equations of motion (2-37) become 



(2-39) 



2 

V = -K.yV - g sin y 

y = K 2 yV + Ve(V) 

y = -j3yV sin y 

r6= V cos y 

Eliminate the independent variable, time, in favor of the variable y by dividing each 
of the other equations through by y. Equations (2-39) then become 

(2-40) 
(2-41) 



and 



dV V 


K lV j g 


dy y 


|3 sin y /3yV 


dy _ 21 - 
dy y 


K 2 e(V) 

(3 sin y 0y sin y 


de_e_ 

dy y 


ctn y 

)3ry 



(2-42) 



These are the equations to be integrated for V, y, and as functions of the vari- 
able y, which in turn is expressible in terms of the altitude h. It should be noted 
that no assumption has been made as to the magnitude of y; that is, no assumption of 
a small flight -path angle has been made. 
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Variation of flight-path angle y .- In order to evaluate the flight -path angle y, 
the following steps are taken. A final, or break-out, condition is assumed at the end 
of the trajectory. This condition could be used to specify, for instance, the altitude 
of the end of the perceptible atmosphere, that is, h b , skipping out, or the end of an 
entry-control mode at some height or flight -path -angle, y., condition which is desired 
by the vehicle guidance. In any event, the subscript b represents a local, final flight 
condition. Then the integral of equation (2-41) becomes 

If , during the interval of evaluation, e(V) has a very weak variation with y, then 

y b K„ 



-cos y 



y ' Y 



Yb e(V), 

— i—i In y 

y P 



y b 



or 



cos y = cos y h + -2(y - y b ) + ^ In L. (2-43) 

In terms of h, equation (2-43) becomes 

cos y = cos y b + -^(e"^ h - e b ) - e(V)(h - hA (2-44) 

For the appropriate values of j3h and j31v 

l(e-*. e -^) = . (h . hb) 

Since K2 » e(V) for most trajectories (5 x 10 -7 as opposed to -1 x 10" 7 as deter- 
mined numerically at the point of maximum dynamic pressure and less above this 
point for a typical Apollo entry), it can be seen that the third term on the right-hand 
side of equation (2-44) can be neglected compared to the second term on the right- 
hand side. 

Variation of the velocity V .- For purposes of evaluating velocity V as a func- 

.n the form 

(2-45) 
(2-46) 
(2-47) 



tion of 


y, equation 


(2-43) is writt 




cos 


y = 


a + 


by 




a = 


cos 


^b 


J b 


and 


b = 


K 2 








Equatio 


n (2-45) 

sin 


ma 

y = 


y be rewritten 




Vi. 


• (a + by) 2 



(2-48) 
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or 



sin r = ||(l - a 2 ) - 2aby - b 2 y 2 



(2-48a) 



Now the integral of equation (2-40) with the proper limits yields 



r v bdv = i r y b/ K i l g 



J v V PJy \siny y 2 y 



dy 



(2-49) 



Again, it may be seen from numerical calculations that the second term on the right- 
hand side of equation (2-49) will contribute very little to the solution because 

L_ » JL by several orders of magnitude. Typically this number is 2.5 X 10" as 

sin y y2 

opposed to 2.5 X 10" 7 for an Apollo spacecraft at the point of maximum dynamic pres- 
sure. Furthermore, the small variations in g and V with y throughout the 
interval of integration will have even less influence on the solution and will therefore 
be neglected. Treating g and V in the second term of the right-hand side of equa- 
tion (2-40) as constant average values (denoted by g and V) and replacing sin y 
with the identity of equation (2 -48a) alters equation (2-49) to 



Ki 



ln^-=~ -i= + 



Vb 



^bdy 



ft' 2 



M k 



(2-50) 



where 



Y = (l - a 2 ) - 2aby - b 2 y 2 (2-50a) 

By the use of integral tables and several trigometric identities, the integral of the 
first term on the right-hand side of equation (2-50) may be put into the form 



ft" 



4 si"" 1 

±b 



a + by 
±1 



(2-51) 



From equations (2-47) and (2-45), 



.-1/-CQS y\_ P 



f^^Asin-iterW 

J ^/f K 2 \ ±1 / 



K, 



i- - cos 



-1/ -CQS y\ 
\ ±1 ) 



(2-52) 



Depending on the choice of signs for the second term of the right-hand side of the 
preceding equation, the integral of dy/VY may become 



f J 



* S-H - y) 



or 



yf K 2^ 2 



dy B (it , \ 

w =± 4( 2+y ) 



(2 -82a) 



(2 -52b) 



In either case, the same answer results when the limits of equation (2-50) are applied 
and when the values of equations (2-52) are substituted into equation (2-50), namely, 
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ln V b ±1^ y W+ 2 ^y b 



0V 

Clearing logarithms and replacing y with e 
following explicit expression for V 

D 



|3h 



and Kj / A^ 2 w® 1 ~J~ yields the 



V 



■L'( r - r b)-^ 



i h -\) 



(2-53) 



Again it should be noted that equation (2-53) is valid only in an interval in which 
the second term of the exponent is always small compared to the first term. This 
second term should be considered an approximate second-order correction to the first 
term by virtue of its derivation. It also follows that, inasmuch as Vjj must always 

be less than V (because of atmospheric drag), the first term of the exponent must 

- /— 2 

always be a positive quantity. Finally, because a value of g/V must be assigned in 

equation (2-50), it should suffice to use g b and V^ if one is solving for V and to 
use the present-state values of g and V if one is solving for V^. It should also 
be noted that this equation is not applicable for the case of zero lift although the zero- 
lift case can be bracketed quite nicely by applying this equation with lift values very 
near to zero. 

Variation of range 8 .- Integration of equation (2-42) for range yields the 
following equations: 



J 



de 



= _ JL r y k ctn y dy 

y 



^rj y 



or 



T^de 



if 1 



a + by 



dy 



0rJ, 



yb^TT 



(a + by) 



2 y 



and 



J 6 ^ \ y fi r J y b # 



At this point, three cases must be considered. First, if a 2 = 1, then 

^y 



,3r\aby/y b 0r 



% 



J- a in-V - C0B r 



Ko 



±1 



= ± 



K 2 r 



1 /sin y 



sin 



rb 



yb 



pr 



(^ - Tb) 



= ± 



K 2 r 



fiKu 



/Y, 



e-~sin y - e "sin y h j ± ±-(y - ^ 



(2 -54a) 
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In most practical cases a 2 is very nearly unity, therefore equation (2-54a) will 
usually suffice. The following possibilities are included, however, for the sake of 
completeness: Second, if a 2 < 1, then 

1Y 



% 



- = ± 



prVl^a 2 



-In, 



JY+ Vl-a^ ab 



iT^i 



p 



^( r " y b) 



Jyb 



0r Vl - a* 



-lnSHLii 



+ Vl - a 2 



ab 



L 



l/TV* 



sin 



+ ln 



^b 



+ /T 



ab 



?b 



^/T> 



± Jrt - n>) 



/3rVl - a 2 



In 



1 - a 2 - aby b + yl - a 2 sin y b / y 



1 - a 2 - aby + yl - a 2 sin y ^ ■ 



1_ 



± i^[y - r b ) 



= ± 



/3ryi - a *- 



»W-*-V 



±U-r h ) (2-54b) 



where z(y,y) = 1 - a 2 - aby + h - a 2 sin y. Third, if a 2 > 1, then 

iy 



% 



= ± 



•ya 2 - 



sin 



If -aby + 1 - a' 

±by 



|3r 



(? - >b) 



y b 



= ± 



j3rVa 2 - 1 



±by / \ ±by b 



b' 



: (y - y b ) 



= ± . 



a / . -1 

'sin 



•£2 



|3(1 - a cos y) 

±Koy 



- sin 



-1 



K 1 - a cos n>) 



±K 2 y b 



h{ y - y b) 



(2 -54c) 

In a more compact form, equations (2 -54a), (2-54b), and (2-54c) may be expressed as 



%- e -h 



r(y,h) ± (y - r b ) 



(2-55) 



where the negative sign on the term in brackets is the proper choice in all cases, and 
where, for the first case when a 2 = 1, 



38 



r(y,h) = ± |-(e^ h sin y - e b sin yj 



for the second case when a fi < 1, 



T(y,h) = ± 



R 



i(h b - h) 



'1 - a - abe 



-0h, 



+ ln 



M 



1 - a' 5 sin 



and, finally, for ?r > 1, 



T(y,h) = + a /sin" 1 



f2TT 



,/3h 



K 2 



(1 - a cos y) 



- sin 



V 



1 - a 2 - abe^ h + \/l -a 2 



sin y , 



j3e 



Kn 



"I 1 



a cos 



(2-56a) 



(2 -56b) 



(2-56c) 



For cases run in negative time (that is, those cases where the sign on time or else the 
sign on all of the velocity terms is reversed), the positive sign in the term on the right 
in equations (2-56b) and (2-56c) is proper. For cases run in positive time, the nega- 
tive sign is proper. This sign-selection convention applies even though time does not 
occur explicitly in the solution in order to keep the down-range angle consistent with 
the flight-path angle and altitude. 



Kepler Formulation of the Entry Equation 

Incorporating Lift and Drag 

Equations (2-2) and (2-3) can be generalized from the foregoing to include the 
effect of lift and drag. This formulation has obvious applications in the theory of 
decay of near-earth satellite orbits where the drag acts as the main perturbing influ- 
ence on a normally elliptic orbit. Consider figure 2-4 and resolve the lift and drag 
forces into their respective radial and tangential components. The radial and tan- 
gential lift and drag coefficients are given by 

C r = C L cos y - C D sin y 

C = -C L sin y - C D cos y 
Equations (2-2) and (2-3) become, in terms of forces, 



(2-57) 



mr - mrB + 



2 m g e r 



e e 



= F(r) 



(2-58) 



mr'e + 2mr0 = F(0) 

From figure 2-4 

F(r) = L cos y - D sin y 
F(0) = -L sin y - D cos y 



(2-59) 

(2-60) 
(2-61) 



39 



L sin 




Local horizontal 



Figure 2-k.- Forces on a reentry vehicle in rectangular coordinates. 
but from equations (2-4) and (2-5) 



cos y = 



sin y = 



re 



(2-62) 
(2-63) 



Hence, replacing L and D with their equivalents, from equations (2-18) and (2-16), 
respectively, 

l.,r2^ c[vQ\ 1 Tr 2~ J'r \_ 1 



F(r) = lpV 2 C L S^j - IpV 2 C D S^j = |pS(r0C L - rC D ) V (2-64) 

which, on substituting for V, becomes 



F(r) = |ps(r6C L - rC D )(r 2 + r 2 6 2 ) 



(2-65) 
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Also, 

Y(8) = - |pV 2 C L s(£) - |pV 2 C D s(f y - ip S (rC L + reC D )v (2-66) 

which, on substituting for V, becomes 

F(6) = - |pS(rC L + r0C D )(r 2 + r 2 f) (2-67) 

Then equations (2-58) and (2-59) become 

.2 m S e r e 2 1 / • • ^(.2 9 -2\ 1/2 

mr - mrfl + ^- = ±pS(reC L - rC D )(r + r 2 j (2-68) 

a/2 



mr 



9 + 2mr0 = - |pS(rC L + r0C D )(r 2 + v 2 § ) (2-69) 



These equations are formally equivalent to equations (2-19) and (2-20). It should 
be noted that, in spite of the familiar form shown here, these equations are coupled 
both in angular momentum and energy because of the lift and drag terms. Of course, 
these equations reduce to the usual Keplerian form if p = 0. Thus the equations of 
motion of a satellite with drag can be compared easily to the equations of motion of 
the same satellite without drag. 

Energy and Momenta Rates 

Although in equations (2-68) and (2-69) the total energy and the angular momen- 
tum are coupled, the same does not apply to the time derivatives of these quantities. 
For example, at any instant of time, the angular momentum is given by 
2 . 
mh = mr 6 (By definition of angular momentum) 

Then the derivative with respect to time is 

• 2- 

mh= mr a + 2mrr0 (2-70) 

This value of mh, however, is the left-hand side of equation (2-69) multiplied by r, 
and hence this equation can be written as 



mh = - ips(c D r 2 + C L rr)(r 2 + r 2 2 ) 



2 -\~u- - ■ -l-v^ ' - " ' (2-71) 

Also, at any instant of time the total energy is 



me r 



2 



1 O &c 

mH = - mV z E_J_ (By definition of total energy) 

or, taking the time derivative, 

. m g e r p 2 . 

mH = mVV + — %f- r (2 -72) 
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Also, the change in energy is equal to the integrated product of the drag and the dis- 
tance traveled. With s as distance and Hq as the initial energy 

m(H - Hq) = - J S |pSC D V 2 ds ( 2 " 73 ) 



but 



so 



ds = V dt 



m(H - Hq) = -J |pSC D V 3 dt -(2-74) 



The energy rate can thus be put in the simple and elegant form 



1 pSCpV (2-75) 



H=_ 2 m 

This is an important conclusion in several respects: First, it is quite general; that 
is, there are no strained or artificial assumptions in its derivation. Second, it 
applies at any point along a reentry trajectory and thus depends only on the density at 
the present orbit. Third, although it applies to reentry equations containing lift, it 
does not, itself, contain any lift terms. And fourth, it is a function of the single 
dynamic parameter V 3 . This is as it should be. It results from the work-energy 
theorem of classical mechanics. 

The corresponding angular momentum equation can be found in an analogous 
manner. By application of the transformations 

v6 = V cos y 

r = V sin y 
in equation (2-71) and cancellation of a few terms, the angular momentum is found to 
be 

h= - iPS m V2 ( C D C0Sy + C L Siny ) (2 " 76) 

It can be seen that the rate of change of angular momentum per unit mass does depend 
on C L as well as on C D . 
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APPENDIX TO CHAPTER 2 

DERIVATION OF EQUATIONS (2-68) AND (2-69) 

The transformations developed in chapter 2 can be handled more elegantly in the 
following manner for those who prefer matrix techniques: 

Equations (2-7) and (2-8) can be written 



m 



cos y -sin y 



Sin y cos y 



Vy 



mf 



— 




1 — — 


sin y cos y 









cos'V 


- m 


g r 2 
& e e 

L r2 J 







(2A-1) 



Denote the first matrix on the left by 
cos y -sin y 



sin y 



cos y 



(2A-2) 



Hence 



A" 1 - 



cos y 



sin y 



-Sin y cos y 



1 



Premultiply both sides of equation (2A-1) by A" to obtain 



(2A-3) 



m 



Vy 



mV £ 
r 



cos y sin y 



-sin y cos y 







sin y cos y 






- m 


COS y 





cos y sin y 



-sin y cos y 



E r « 
6 e e 



(2 A -4) 



Thus 



r 2 



mV / 2 2 ^ m £p r p 
mV = — — (^-sin y cos y + sin y cos y) ^ 



sm y 



which reduces to 



mg r * 
mV + — £- sin y = 



(2-11) 



Also, 



mV 2/ 2 3 N m g p r „ 
mVy = — — -\cos y sin y + cos yy £- . ■ cos y 
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APPENDIX TO CHAPTER 2 - Continued 



which reduces to 



. mV 2 m g e r e 2 

mVy = HLL- cos y - e „ e cos y 



or with 



m g e r e 

mg = __ 



£4 



mVy - mgj- - 1 cos y = 



(2-15) 



From equations (2-19) and (2-20), 



m 



V 


= -mg 


sin y 


- ipSV 2 
2 


C D 


v(y- V 7 y ) 




cos y 




_" C L_ 



(2A-5) 



Premultiply both sides of this equation by A 



m 



cos y -sin y 



sin y cos y 



V 



„/• V cos y 



= -mg 



cos y -sxn y 



sin y cos y 



sin y 



cos y 



ipSV 2 
2 



cos y -sin y 



sin y cos y 



'D 



(2A-6) 
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APPENDIX TO CHAPTER 2 - Continued 



which becomes 



m 



V cos y - v(y - V °° S y W n y 



V sin y + V 



( f _vcosy) 



cos y 



= -mg 



Sin y cos y - sin y cos y 



2 2 

sin y + cos y 



|pSV 2 



Cp cos y + C L sin y 



Cq sin y - Ct cos y 



Equations (2-4) and (2-5) are 

r6 = V cos y 

r = V sin y 
Taking the time derivatives of these equations 

re + r'e = V cos y - V sin yy 

r = V sin y + V cos yy 
Substituting these values back into equation (2A-7) 



m 



re 



3 + T6 + v cos y sin Y 



2 2 

i? - v cos y 

r 



= -mg 







- |psv 2 



Cj3 cos y + Cl sin y 
Cj) sin y - Cl cos y 



(2A-7) 



(2-4) 
(2-5) 

(2A-8) 
(2 A -9) 



(2A-10) 



Then, substituting various combinations of equations (2-4) and (2-5) to eliminate the 
remaining terms in V and y this equation becomes 



°ni¥W?, 



m 



re + r'e + r§ 











= -mg 




-IpS 
2 


r - re 2 




1 




1— . _J 




- — ' 





D\? 



c d1If)- c l(y 
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APPENDIX TO CHAPTER 2 - Concluded 

which is equivalent to equations (2-68) and (2-69) 

/ vl/2 

mr - mr0 2 + mg = IpS (c L r6 - C D r)(r 2 + r 2 6 2 ) (2-68) 



1 fn 

p2fl + 2mrre = - Ips(c D r 2 6 + C L rr)(r 2 + r 2 2 ) (2-69) 
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CHAPTER 3 
RELATIVE MOTION IN RECTANGULAR COORDINATES 
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L-73-3011 
The third stage of the Saturn V used to launch Apollo 9 just prior to docking 
for removal of the Landing Module. This photograph was taken in earth 
orbit. 
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CHAPTER 3 

RELATIVE MOTION IN RECTANGULAR COORDINATES 

An important category of orbital motion is relative motion between two bodies 
which move about the same gravitationally attracting mass; for instance, it may be 
desired to rendezvous one spacecraft with another or with a permanent space station. 
Alternatively, it may be desired to keep track of a spacecraft by knowing its position 
at all times%ith respect to some simpler, better known orbit such as, for example, a 
circular orbit of the same period. Precise and accurate knowledge of relative dis- 
tances and speeds can lead to considerable savings in terms of payload into orbit for 
a given size launch vehicle and, therefore, is an important cost item in modern space 
travel. The need for accuracy has meant that actual orbit determinations and rendez- 
vous procedures are presently carried out numerically by means of large digital com- 
puters and ground tracking arrays. An understanding of the mechanics involved, how- 
ever, that is, an understanding of how spacecraft can be expected to behave in a given 
situation in any given reference frame, remains an indispensable part of the design 
process. The purpose of this chapter is to convey an understanding of the way space- 
craft appear to move when viewed from some other vantage point than the center of 
the earth. 

In this chapter for the first time a fixed inertial reference frame is not used and 
instead a system is considered in which the coordinates are, themselves, in motion 
with respect to an inertial frame. The material in this chapter is, in fact, developed 
from the classical theory of moving coordinates and in every case is merely a special 
application of the theory. Finally, it is remarked that relative motion is the simplest 
and most tractable case of three-body motion, that is, three bodies, only one of which 
is gravitational, and the only case that can be handled adequately in a practical way 
without extensive numerical calculation at the present time. For purposes of this 
book, this is the nearest to the complicated problem of three -body motion that will be 
considered. 

Rectangular Coordinates 

Rectangular coordinates are covered first for three reasons. The first is his- 
torical, since the first rendezvous equations were presented in rectangular form. 
The second is pedogogical. For the beginner accustomed to thinking in rectangular 
coordinates, these are usually the easiest to understand. Third, rectangular coordi- 
nates are as accurate as any other system provided the bodies under consideration 
are in close proximity. Thus, this approach is still advantageous for some problems 
in rigid -body dynamics, and, in fact, rectangular coordinates will be used in the fifth 
chapter of this book when rigid bodies are discussed. Rectangular coordinates are 
not superfluous at the present point in development either, however, since shell coor- 
dinates develop rather naturally from this form. It will be seen, in fact, when the 
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Clohessy -Wiltshire equations are developed in this chapter in rectangular coordinates 
and in the next chapter in shell coordinates, that the equations are the same to first 
order in both systems. 

The coordinate system in this approach is shown in figure 3-1. The system 
consists of rectangular coordinates centered on one of the bodies in orbit. In the case 
of rendezvous this body is usually the target vehicle. 



-Target vehicle 
R 




Rendezvous 
vehicle 



Figure 3-1.- Rectangular-coordinate system centered on a body in orbit. 

In equation (1-9) we gave the general expression for the Lagrangian in a non- 
inertial frame of reference without translation. The equation was immediately 
reduced to the form of equation (1-10) by the assumption that Q, is zero. We now go 
back and pick up the argument assuming that fi is not equal to zero. Equation (1-9) 
is reproduced here for the sake of convenience 

-2 - - - 1 trr -\ 2 



L = - mV 2 + mV • fi X r + i m(o x r) 
2 2 ' 



r 



In the case of moving coordinates (fig. 3-1) it is assumed that 
5 = k0 



(1-9) 



(3-1) 
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and that 

r = Tx + T(y + r s) + £ z 
Hence , it is found that 

fix?= -Te(y + r s ) + i 6x 

Also, for rectangular coordinates, the relative velocity is 

t v = Ti + T(y + r s ) + ^ 

Then 



(3-2) 
(3-3) 

(3-4) 



V • fix r = 



x(y + r s ) -x(y + r s ) 



and 



(3-5) 



(H x r ) = (y + r s ) 2 + x 2 



(3-6) 



In rectangular coordinates J^-. is given by 

r 






-% r e 



~o 2 21 1 / 2 

x + (y + r g \ + z 



(3-7) 



Hence, the Lagrangian becomes 



L = -z m 



.2 /. . N 2 .2 

x + (y + r s ) + z 



+ ra< 



x(y + r s ) -x(y + r s )" 



1 a2 

+ ±m0 

2 



(y + r B ) 



2 + x 2 



mg e r^ 



x 2 + 



(y + *sf + z 2 



1/2 



(3-8) 



Taking the appropriate derivatives, the exact differential equations of motion are 
arrived at in rectangular coordinates 



mx - m(y + rSe - 2m(y + r s )e - mx0 + — 



mg„r„ "x 
b e e 



-,3/2 



= 



(3-9) 



•2 m S e r e (y + r s) 
my + mx'e + 2mx0 + mr s - m(y + r s \6 + - 372 = ° 



x" + 



(y + r s y 



+ z' 



(3-10) 



mz + 



mg e r e z 



x 2 + 



(y + r s) 2 + z *\ 



3/2 



= 



(3-11) 



51 



If there are nonconservative forces such as thrust or drag, terms incorporating these 
forces may be added on the right-hand side of these equations. It will be assumed for 
the present, however, that there are no nonconservative forces acting on the vehicles. 
By factoring out the term r s 3 in the bracket, it can be established that 



1 



x 2 + | 



\£j , yCl T* "5 



(y + r s) 

Then, equations (3-9) through (3-11) become 



^£) 2 *M 



r s* r s 



3/2 



(3-12) 



.2 m Se r e 

mx - m/y + r s \0 - 2m(y + r s )0 - mx0 + ^ — 



i + 2_f + *L + -4 



S/ t ^ r ^ 

a/ r s r s 



-3/2 



(3-13) 



.2 m g e r e 2 ( y+r s) 

my + mx6 + 2mx6 + mr s - m/y + r g \0 + ^ 



1 I_\ 2 + EL + Zl— 
+ r s) r g 2 r s 2 



1-3/2 



mg p r " 
mz + Z-Z-z 



1 XJ\ 2 + 2EL + 3— 
r s) r s 2 r g 2 



-3/2 



= 



= 
(3-14) 

(3-15) 



This result may be specialized by assuming that the origin of the coordinate system is 
in circular orbit about the attracting planet. Then, denoting 6 by the constant 
angular rate u> and observing from equation (1-54) that 



!iV = w 2 



(3-16) 



Also, if 



r s = 



r s = 



and also, if 



e = o 



Then 



o 2 

mx - 2mwy - morx + mco x 



!+:£- + 



r s y r 2 r 2 



-3/2 



* 



sj 



(3-17) 



my + 2mwx - mw 2 /y + r g ) + mw Z (y + r s \ 



1 + xx + ^ + z ' 



V7? 



-3/2 



(3-18) 
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mz + mu z 



1 + 



y 



-3/2 



= 



(3-19) 



Numerical solutions to equations (3-17) and (3-18) are presented in figure 3-2 
which is modified from reference 5 for the special case where the z coordinate has 
been held zero. In this figure x is plotted against y for a vehicle which is ejected 
from the origin at various ejection angles but at a constant speed of 3 meters per sec- 
ond. The dashed lines show the locus of positions of the ejected vehicle for various 
angles of travel about the center of the earth. The effect of trajectory curvature is 
quite evident from the spiral nature of the resultant figure. 

8 

V Q = 3 m/sec 



4- 



y, km 



-6- 




75/ 

/ 9, deg about 
/ center of earth 



-8 



-6 



-2 




x, km 



Figure 3-2.- Trajectories as measured in a coordinate system above the 
surface of the earth of a number of point masses (which at t = 
were at the center of the coordinate system) with different initial 
velocity components relative to the moving coordinate system. 
Dashed contours denote the motion of the masses at subsequent posi- 
tions of the origin of coordinates. 

The last terms on the left in equations (3-17), (3-18), and (3-19) suggest that an 



expansion of 



\ + r 



N*IH) 



fl-3/2 



in a power series might be useful in 
canceling terms. Such an expansion is facilitated by the fact that under most practi- 
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cal circumstances x, y, and z are all much smaller than r s so that reasonably 
rapid convergence can be expected. Expanding (l + £) and defining the small 
term v as 



The term to be expanded becomes 



(3-20) 



-3/2 



Then, keeping three terms, an expansion of the form 



/* \-3/2 

(1 + v) 



(l + ,)- 3 / 2 «l-|, + ^v 2 + ... 



gives 



' 1 + $ + 



ft) 2 + (ft) 



-3/2 



-1-1 



® 2+ ft) 2 + ft) 



15 
8 



2^- + 



;^) 2+ (^ 



(3-21) 



By squaring the last term on the right and retaining terms correct to second order, 
equation (3-21) becomes 

-3/2 



(!♦#♦$♦(£ 



1 3y_ + 6yf_3 



s r, 



2 2 



(¥ + ® 



(3-22) 



Making the appropriate substitutions in equations (3-17), (3-18), and (3-19) and can- 
celing terms of higher order than the second gives 



mx - 2mwy - 3mw 21 = 
*s 



• „ 2 o « V 3 m u) 2 x 2 3 m cA 2 . 
my + 2mu>x - 3mw y + 3m -£_ - 5 m — - g- 2 m r 2 



= 



mz + mu z 



2~ _ 3mw 2 ^=0 



(3-23) 



(3-24) 



(3-25) 



Clohessy -Wiltshire Equations f 

Simplification may be continued by neglecting the squared and cross-product 

terms since these terms are small. Then 

(3-26) 



mx - 2mcuy= 



my 



+ 2mo)x - 3mw y = 



(3-27) 
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mz + ma; z = (3-28) 

These, then, are the first-order relative -motion differential equations. 

A first integral of equation (3-26) can be obtained directly since both terms are 
time derivatives. Then x from the solution of this equation can be used in the sec- 
ond term of equation (3-2?) to obtain a complete solution in terms of y. Finally, this 
solution can be applied in the first integral of equation (3-26) to solve this equation 
completely. Equation (3-28) is, of course, just simple harmonic motion and so can be 
solved directly. Thus, these three equations have the solutions: 



2 2 



^o 

0) 



Sj J sin wt - 2 -£ cos wt + I 6y Q - 3 -g-Jujt + 2 






(3-29) 



2^> 
as 



3yJcos«t+^ S inwt + 4y -2§ 



(3-30) 



z„ cos wt + -r~ sin wt 

o w 



(3-31) 



These equations give the position of an arbitrary vehicle relative to an origin which is 
itself orbiting in a circular orbit. If, for example, a target satellite is located at the 
origin, then x ,y ,z and x ,y ,z could be interpreted as the initial position and 
velocity, respectively, of a rendezvous vehicle. Equations (3-29), (3-30), and (3-31) 
then are first-order expressions for the subsequent motion of the rendezvous vehicle 
with respect to the target. It follows that, if it is desired to achieve interception at 
some specified time Tj, it is merely necessary to set x, y, and z equal to zero 
and solve for the present velocity needed to achieve these conditions. The instanta- 
neous relative velocity components required are given by the following expressions: 



to 



3 (sin WTjj + Yo^wtj sin wtj - 14 (l - cos wtS 
3o)Tj sin WTj - 8(1 - cos wtA 



(3-32) 



y c __ 2x (l - cos wtj) + y Q (4 sin WTj - 3wTj cos wTj) 



w 



3wt- sin (*>t* 



8/1 - COS WT, 



(3-33) 






-2, 



tan wt. 



(3-34) 



These are the rendezvous guidance equations first published by W. H. Clohessy 
and R. S. Wiltshire in 1960 (ref. 6). They are. frequently referred to as the Clohessy- 
Wiltshire equations, and their influence on the development of modern analytical astro - 
dynamics has been substantial. It should be mentioned that these equations were first 
used to demonstrate the technical feasibility of earth orbital rendezvous by Eggleston 
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and others (refs. 5 and 7). The equations were actually incorporated in -the Gemini 
spacecraft rendezvous computer (although in shell-coordinate form) for the first suc- 
cessful attempts at earth orbital rendezvous. Subsequent hardware developments 
have tended to rely more heavily on numerical integration with digital computers due 
to the high accuracy attainable, but these equations are still useful for analytical and 
developmental work. 

The great drawback of the Clohessy -Wiltshire equations, as might almost seem 
obvious from this derivation, is that they break down quickly. The exact differential 
equations describing relative motion, equations (3-9), (3-10), and (3-11), are nonlinear 
and cannot be solved analytically. A number of approximations were made to yield 
equations (3-26), (3-27), and (3-28), which are linear, and the solutions are given by 
equations (3-29), (3-30), and (3-31), which in turn lead to the Clohessy -Wiltshire equa- 
tions, (3-32), (3-33), and (3-34). Because of the simplifications made, the Clohessy - 
Wiltshire equations are accurate only for relatively short time periods and for dis- 
tance not too remote from the origin. The nonlinear effects, which were ignored in 
order to obtain analytical solutions and which, though small near the origin and for 
short time periods, become larger and influence the solutions more 'at greater dis- 
tances and longer times. This loss of accuracy with distance and time can be seen in 
figure 3-3 where the solutions of equations (3-29) and (3-30) are compared with exact 
numerical solutions done on a computer. 
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Figure- 3-3- ~ 'Loss' of accuracy of the first-order relative-motion equations 
/^t'i <i^ J ii- r if-. ' -i >- v rtii t-lar - -mi." 1 - -~ r • ■> _iina,te sys- 

tem in circular orbit above the earth. 
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Yarious attempts have been made to improve the accuracy of both equa- 
tions (3-29), (3-30), and (3-31) and, also, (3-32), (3-33), and (3-34). In the next por- 
tion of this chapter two subsequent attempts at improving the relative motion equa- 
tions (3-29), (3-30), and (3-31) will be discussed. The distinguishing feature of these 
two improvements is that they are in rectangular coordinates. A similar attempt in 
nonrectangular coordinates will be taken up in chapter 4. 

London's Improvement of the Relative -Motion Equations 
A sefeond-order approximation to the solution of equations (3-17), (3-18), and 
(3-19) can be obtained (ref. 8) by referring back to equations (3-23), (3-24), and (3-25). 
Note that these equations contain terms in the gravity field expansion, equation (3-22), 
which were neglected in the previous linear development. 
Let 

X = Xj + x„ 

y = j 1 + y 2 

z = z 1+ z 2 

where the subscript 1 denotes the first-order solution and the subscript 2 denotes 
a second-order correction. Substitute these values back into equations (3-23) (3-24) 
and (3-25) 



mxj + m£ 2 - 2mojj l - 2maiy„ - 3mw 2 Li JLlLJ__21 = o 



my j + my 2 + 2u)k t + 2wx 2 - 3ma> 2 j 1 - 3mw 2 y, + 3mw 2 i!i_LM_ 



(3-35) 



-fm^fcL^-Wkp)!.,, (3 _ 36) 

r s s 

mij + mig + mw Zj + raw z, - 3ma) 2 i-JL. — Mi_J____2J _ Q (3.37) 

s 

Immediately it is seen that, since x 2 , y 2 , and z 2 are to be thought of as small 

correction terms, higher-order terms involving x 2 , y 2 , and z 2 may be ignored. 

Hence, in equation (3-35) 

2 
3moj , 9 x 1 y 1 

- -y— (xjyj + x 2Fl + Xl y 2 + x 2 y 2 ) * -3m w^ -±± ( 3 _ 38) 

s 
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3m cu 2 / 2 , -„ „ , „ 2\« w.^Ii 



M y i 2 + ^ + y 2 2 ) s 3ma)2 t~ (3 ~ 39) 



Likewise, 



3mc- 2 ' » "^ ' '" 1 

2r s 



-( X l 2 + 2 V2-2 2 )^|- 2 i (3 " 40) 



■ z 2 

32^1(z 2 + 2z 1 z 9+ z 9 2)«-i mw 2 1 (3-41) 

2r . V 1 j. <2 ^ / 2 '-s 



Finally, 



Then equations (3-35), (3-36), and ,(3-37) become 

9 Xi ¥ 1 
mL + rain - 2mcoy, - 2mujy 2 - 3mo> — — = . (3-43) 



s 



my. + my„ + 2m«x, + 2mwx 2 - 3mw y^ - 3mw y g + 3mw -^ 






x 2 z 2 

. 1 mo) 2 J_ - i mJ J- = (3-44) 

2 r s 2 r s 

2 
:«» ■ *" 2 2 3mti? t% /o ji£\ 

biz, + mzg + mi) z- + mw z 2 - — j — yjZj = , (d-4o) 

■S 

Subtracting the first-order equations (3-26), (3-27), and (3-28) from these, respec- 
tively, yields 

mX 2 - 2m W y 2 = *f^- x^ ' (3-46) 



my 2 + 2mo;x 2 - 3m W 2 y 2 = - ^L^ 2 - f (x^ + z^j 



(3-47) 



.2 



mz 2 + m^ 2 z 2 = ^L y lZl (3-48) 

Thus a set of linear differential equations with constant coefficients has been obtained 
for the second-order corrections. The right-hand sides are known functions of time 
since each term can be computed from equations (3-29), (3-30), and (3-31). It is sig- 
nificant that the out -of -plane and in-plane terms are now coupled. We shall find that 
this coupling continues when we investigate other types of coordinate systems later on. 
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Equation (3-46) becomes 

2 

"s 



mx„ - 2mwy„ = — — — (A« cos o)t + B-, sin cut + C« cos wt sin tut 



+ D, wt cos u)t + Ei wt sin wt + F.. cos o>t 



+ Gj sin wt + Ij wt + J A (3-49) 

where 



~4i o y Q 6y Q y o 



4i y Q 6y y 
B __ _ -—57— 
0)* 



V. M ^. + 1 ^> 2 ^ 



to 2 w ° w 2 



Vo 6 x 2 o 

D ,21 ^ . _o_ . 18y 2 



6y y Q 3i Q y o 
El -_ — 



Fl = 8 ^2-14^£ + 2^2^x y 



. 2 • • -2 

r - 9 F ° 4. ?R X ° y ° + X ° y ° 24v 2 ft X ° 
Gj - 2 —j- + 28 -—- + -_- - 24y Q - 8 -g- 



Il = 2 4y o 2_ 2 4 X ^2 + 6 X ^ 



2x o x o , ,Vo a Vo 

J i = - -w- + 4x y Q - 4 -j- + 8 — 
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Equation (3-47) becomes 



"[ A 2 



my. + 2muL - 3mco y„ = — — |A„ cos cot + B„ sin tot + C„ cos tot sin tot 



where 
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+ D 2 ^t cos wt + E2 wt sin wt + F9 cos wt 



+ G2 sin wt + Ho(wt) + In cot + J2 



(3-50) 



A 2 



4x 



co' 



2y 2 z 2 
12 _lJi + 9y 2 . _fo o_ 



k 2 ,0^0 



CO •'o 



• 2 "2 



18y 



2 ^o 
~2co 2 



Vo y o y o z o z 



09 = 12-^-18 M 

* , ,2 co co 

CO 



_o_o , 1<5 y o y o 

CO 



D 9 =-6-^ +12 ,, 

& _£ it) 



E 2 = 2 \ 2 

\ CO" 2 



6x ~ x o y o 



21^_H + 18y o 



* 2 • 2 
F - 28 ^2l° - 24v 2 R X ° + 4 Y ° 4- 9 X ° Y ° 



G , o -12 f<^ - 4 ^ + 20 ^ + 6x,v 



CO 



CO 



CO o'o 



„ _ 9 x o 18 X o y o 1R 2 
H 2 _-_ + l8_-18y o 



x o y o . o x o x o , o y o y o 



I 2 = 6^ + 3 



CO* 



CO 



12 



CO CO 



^rJr 



w co 2 co2 2 



J 9 -16y'- 16-^ + 4-^--^- -2 



Equation (3-48) becomes 



2 3mw 2 /. 2 l4 . „ „. 2 



where 



1W / 

mz, + moJ^Zg = ^~ — ~( A 3 cos"wt + B„ sin^wt + Cg cos wt sin wt 

S 



+ Fg cos wt + Go sin art] (3-51) 



A, = 2 -2-2. - 3y z 

o 0) j O C 



¥ Z 
B = ° ° 

w 



X Z V Z V z 

n o oOO^OO 

C 3 - 2 -j- - 3 -^- + -j5- 



F,=4yz - 2 ° ° 



3 y o o u> 



G 3 = 4^2-2^ 



Solutions of the form 

x, = a + a* wt + au sin w t + a, cos art + a, sin 2 art + a,- cos 2 art 

+ agart sin wt + a„wt cos art (3-52) 

2 

y 2 = |3 + ^jWt + j3 2 ( wt ) + J 3 3 sin wt + ^4 cos wt + % sin 2a> t 

+ (3g cos 2 art + j3 ? wt sin wt + ;3gWt cos wt (3-53) 

Zg = y Q + y, sin wt + y 2 cos wt + Y% sin 2 art + y, cos 2art 

+ ygO>t sin wt + ygWt cos wt (3-54) 

are sought. 

The reason for this particular form is evident from an examination of equa- 
tions (3-49), (3-50), and (3-51). If these equations are converted to a multiple -angle 
formulation, it can be seen that compatibility upon differentiation demands just this 
form so that a balance of terms between the right-hand side and the left-hand side is 
possible with no degeneracy. The constants a n , /3 , and y are functions of the 
initial values x ,y ,z o and x ,y o ,z o , of course. The method of solution which will 
be adopted here is to substitute these solutions along with the appropriate time deriva- 

61 



tives back into equations (3-49), (3-50), and (3-51) and then to compare coefficients 
term by term on the right-hand side. The resulting set of simultaneous equations will 
then yield the appropriate values of a^, fi^, y i? and so forth, in terms of A|, B|, 
and Cj. 

The first derivatives of equations (3-52), (3-53), and (3-54) are 

X9 

— = a. + /a, + a„jcos wt + (a* - a^sin c*jt + 2a a cos 2 art 

- 2df, sin 2wt + auwt cos art - a„wt sin wt (3-55) 

^2 



w 



z 2 



/Sj + 2p 2 tat + fp„ + (3q)cos wt + (Prj - |3 4 )sin wt + 2j3g cos 2art 

- 2j3g sin 2 art + |3„wt cos wt - |3gWt sin wt (3-56) 

"aT = f r l + y 6) C0S Wt + f r 5 " y 2) sin Wt + 2y 3 C0S ^ wt " 2 ^4 sin ^ wt 

+ y g wt cos wt - y g wt sin wt (3-57) 

and the second derivatives are 

— - = -(a 2 + 2aJ)sin wt + (2a„ - o<q)cos cut - 4a^ sin 2ort 

- 4« 5 cos 2art - a<gWt sin wt - a„wt cos wt (3-58) 

y 9 

— = 2^ 2 - /j3 3 + 2/3g\sin wt + ( 2/3 ? - ^cos wt - 4j3 g sin 2wt 

- 4^g cos 2wt - j3-wt sin wt - /3gWt cos wt (3-59) 

— = -/y. + 2yg\sin wt + (2y 5 - yg^cos wt - 4y„ sin 2wt 

- 4y^ cos 2wt - y 5 wt sin wt - y„wt cos wt (3-60) 

When these substitutions are made in equations (3-49), (3-50), and (3-51) and coeffi- 
cients of like terms are equated, the following sets of simultaneous equations result; 
From the x-equation 

3 p^lJ = 8(a 5+% ) <3-81b) 
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3C 

i=8(0 6 -a 4 ) (3-61c) 



r s 



3D 

_1= -(a 7 + 2|3 7 ) (3 -6 Id) 



s 
3E 



r s 



1= 2/3g - a 6 (3-61e) 



3F 

"s 

3G 



^-1=2^ -a 3 -20 3 -2/3g (3-6 If) 

1= ~{a 2 +2a 1 + 2p 1 - 2j3 4 ) (3-61g) 



r s 



31, 

^=-4/32 (3-61h) 



s 



From the y -equation 



o /Art -i- Box ojrt 

|Mr— ^) + T^ = 3 ^o " 2a l + 2 h (3 " 62a > 

K~Sr7 = 7 ^ 6 ' 4a4 " 3|35 (3 " 62b) 

3C 

T - 2 -=2(4a 5 + 7 J 3 5 ) (3-62c) 

s 

3D 9 

^=2(2^ 8 + a 6 ) (3-62d) 

s 

3E 9 

ir ^=2(a ? + 2|3 7 ) (3-62e) 

s 

3F 



T ±=2(2^-fi 1 -a 2 -a 1 } (3-62f) 



3G 9 

■— = 2(2^3 + |3 8 + a 3 + a 6 ) (3-62g) 

s 

3H ? 

^-1=3^2 (3-62h) 



3I 9 

-4=3^ (3-62i) 
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From the foregoing equations the independent solutions are obtained directly: 



l/ A 2- B 2\ 7C1 



a 4 = 2l r £ 



8 r s 



(3 -63a) 



a c = 



7/ A l " B l\ 1 C 2 



5 " 8 V r .= 



2 r £ 



(3 -63b) 



a 6- 3 lr7 + r s + r sy 



(3-63c) 



„ _ J F 2 1 E 2 2G 1 
a 7- 3 lr7"^^- 



2 r s r s 



(3 -63d) 



*1 " 2 r s 



(3-63e) 



3^ 

4 r 



H, 



s s 



(3-63f) 



% = i 



1 C 2 lf A l- B ll 



(3-63g) 



l/ A 2- B 2^ 



Jg __, 



2 r « 



(3-63h) 



h 



3^2 3G 1 3 F 2 
2 r s r g 2 r s 



, Gp 3F| 3E< 
__ + __ + _ 

^ s s s 



(3-63i) 
(3-63J) 



In addition to the foregoing linearly independent solutions, the following linearly 
dependent sets of equations are obtained: 



20, + a- = 



3G 3F 



1 



3F 2 3GJ 

2 h~ a 2 = ~ s — r^~ 



(3-64) 
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From the z -equation, 



3 A 3 + B 3 



Y o 2 ' r c 



1 C 3 



J4 



^5 



^6 



2 r c 



1 A 3 " B 3 



(3-65) 



3^3 
2 r s 



l2l 

2 r s 



The terms, y^ and y^, are self -canceling due to linear dependency and hence 
cannot be determined at this point. To get y^ and y 2 one nas to aPPty the initial 
conditions that at t = 0, z, = 0, and z„ = 0. Then 



ri = -(2y 3 +y 6 ) 

y 2 = -(y + ^) 
Integrating the x-equation, equation (3-49), 

3mw fesin 2 «t + 



(3-66) 



mxg - 2mo>y2 



-O.* ~ ID-* 

sin wt cos cot - Ej cot cos cut + D, cot sin cot 



r s \2 
+ (D.^ - G^cos wt + (E^ + F^sin cot + J-^ + i(A. + BA 



cot 



+ —(cut) H-Kj, 



(3-67) 



where K-. is a constant of integration. 

It can be seen very easily by inserting the conditions x 2 = y 2 = at t = that 

(3-68) 



or 



D x - G x + K x = 



K l = G l " D l 



(3-69) 
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From the integration of the x-equation the following simultaneous set of equations is 
obtained: 

|r7= 2 (%- a 4) 



3f A l-V 



'i a 5 + h) 



3E. 
^s" 

3D 



= 2 h - a 6 



3( Dl - G,) 



= a 2 + « 7 - 2(3 4 



3 ( E 1 +F l) 



3 /2J 1 + A x + BjN 



3I 1 

-J- = -2 h 



= -20. 



3K x 

's 1 s 



r. 1 '^ -■>.♦? 



«1 " 2 o 



(3-70) 



From the last of equations (3-70) and the first identity of the y -equation, equa- 
tion (3-62), a solvable simultaneous set is obtained, since 2 is known already: 



o /An + Bn\ 3J 



(3-71) 



r7t G l- D l + T) =a l- 2/3 o 



(3-72) 



Hence, replacing 2 b y H 2/ r s' and eliminating a^ 

3 /A 2+ B 2 \ 3J 2 2H 2 6 Gl 6 Dl 3 Cj 



^o ~ " 2\ r, 



r s r s r s r s 2 r s 



(3-73) 
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Then 



"i--fC 



A 2 + B 2 



, 6J 2 4H 2 9 |G D + . 



(3-74) 



Once /3 is known, one can apply initial conditions at t = to get, from y 2 = 
at t = 0, 



^4 = "Co + %) 



(3-75) 



Hence 



# 



3G, 3F„ 
, 2 _ 



a n = -^-- F ^+2^ 4 



s ~s 
and from y„ = at t = 



?2 



3 =-(01 + 20 5 + 8 ) 



(3-76) 



(3-77) 



Hence 



3G 9 3F« 

a 3 = Tf + Tr- 2 03 

o S 



(3-78) 



and, finally, 



o=-{ a 3 + a 5) 



<x„ = 



(3-79) 



When these equations are solved simultaneously (which is a very tedious and time- 
consuming process), the resultant values are 



«o = 3(x y -^ + fy ^ + I Zo !2 



(3-80a) 



2x 



o 2 i y o . iv „¥o , x y 



„ - * U 2+liv 2 + I z 2 , Z2_ 1 lo_ + i f£_ _ 
ff l - 3 \ X o + 2 y o + 2 Z o + w 2 + 2 w 2 2 w 2 



CO 0) 



(3 -80b) 



« 2 = 36 7vV- 10 x 3x ^" 2 l" 3x ° 2 " Zo2 " 2 7 {3 " 80c 



2 v 



3=-Vo + -J- 6 ^- 2z o 



"o 



(3-80d) 



1 [ l y o 2 x o ^V, 9 2 1 Z o 1 2) 

«4 = 44^-^ + 3 - y o-i y o ngr-iV 



(3-80e) 



l 5"r: 



i!2 3 2 !o + i fo 

2 w \ o w J 2°o) 



(3-80f) 
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a G = 



-4k 4 



(3-80g) 



- 3 l 2 y -lA 2 ^- 3 y 



(3-80h) 



68 



/3 = 3 



**X.i,*-&-%..L*X' 



w 



2 ' 2 Y o " 2 w 2 " w Y o ' 4\"° 



0)' 



(3-81a) 



^l = -^ X o + 2 ^A 2y o-3 



(3 -8 lb) 



^2 



-±-l 2y -22 



(3-81c) 



?3 = 12 yo¥ +6 Vo- 7 ^t- 3x "- + z "- 

O) 



2 + z 22 

0) 0) 



(3-81d) 



£ 4 = - § x Q 2 



5i o 2 1C .. 2, 2y o 2 , 18 Vo *<,* 12 



O) 



2 15 V + ,,2 



W 



O) 



O)' 



2 2° 



(3-81e) 



y /2x„ \ z 2; 

o o „ I o o 



5 o)\o) y o/ 2w 



(3 -8 If) 



p 6 " o y o + , .2 "2 ,..2 4 , .2 " 4 ° w 



cu 2 2 ^2 4 w 2 4 



(3-81g) 



3 7=- 3 (~-%)( 2 yo-^ 



(3-81h) 



*-»£k-£ 



J 8 



(3-811) 



r o ~ 2 



y z / x 

12.2 + z (2-°-3y 

O) o\ w 'o 



(3 -82a) 



V Z X Z V 

v - 3 L2_2 - 12-2 + z 12 
1 w w 2 o o) 



(3 -82b) 



^2 



y,z z x„ 

- 2 -j- - 2 -is- + 3 y z o 



y 3 =- 



z„/ x„ 

— 2-2-3y n 



+ z. 



w 



(3 -82c) 



(3-82d) 



r 4 = "? 



w 



-3y, 



y- = 3z I2y„ - ^ 



W; 



z y 



(3-82e) 



(3-82f) 



^6 




(3-82g) 



The successive -approximation procedure employed here could, of course, in 
principle be carried to any degree of accuracy desired. It should, however, be clear 
from all the effort needed to extract an acceptable set of solutions to even this accu- 
racy that the payoff is hardly worth it beyond this point. There is an advantage still to 
be gained, however, by assuming noncircular orbits. This procedure will be followed 
in the next section so that the results of these two approaches can be compared. 



Analysis of Anthony and Sasaki 

London's perturbation solution can be generalized (ref. 9) to include the situation 
where the origin of coordinates is, itself, in an elliptic orbit. The motion of the origin 
is approximated by means of power series in the eccentricity of the orbit. This solu- 
tion is obtained in terms of the independent variable, time, by a method of differential 
corrections accounting for terms of second order in relative distance and first order 
in the orbital eccentricity of the moving origin. 

As a starting point, rewrite equations (3-9), (3-10), and (3-11) in the following 
form, where some terms have been transposed and where the term -mg r 2 /r 2 has 
been subtracted from both sides of equation (3-10). Thus 



.0 

mx - my0 - 2my0 - mx9 + 


mg e r e 2 x 




(3-83) 


x 2 + (y + r g ) + z 2 


3/2 ~ mi " s + 2mr s e 


.0 

my + mx0 + 2mx9 - m.9 y + • 


m g e r e 2 ( y + r s) 
x 2 + (y + r s ) 2 + z 2 


9 

mg r * 
& e e 

3/2 r 2 

s 




„ mg r ^ 

= -mr„ + mrJ H_E_ 

S S r 2 


(3-84) 

1 
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l 



mz + • 



mg r 'z 



o 9 9 

x* + /y + r g r + z* 



3/2 



= 



(3-U) 



The right-hand sides of equations (3-83) and (3-84), however, are seen to be the same 
as the left-hand sides of equations (2-3) and (2-2), respectively, except for sign and 
are thus equal to zero. Then 



mx - my0 - 2my0 - mxd + 



• " 'A 

my + mx9 + 2mx0 - m0 y + 




= 



(3-85) 



mg r l 
& e e 



= 



(3-86) 



9 9 9 

x + (y + r g \ + z* 



Thus, by the use of an ingenious trick the terms in r and r have been replaced 
by an extra potential term. Also, 



mz + 



mg r 6 z 
& e e 



x* + 



(y + r g )2 + 



3/2 



•= 



(3-11) 



It is convenient to nondimensionalize in terms of the semimajor axis a by 

x = 5 



?-J 



and 



z = l 



P = 



(3-87) 



or 



■-® 3/ wT> 



T = p 3/2 Wt 

In terms of these variables, equations (3-85), (3-86), and (3-11) become 



X - Y0 - 2Y0 - 2 X + 



x 2 + (y + p) 2 + z 2 ] 



3/2 



= 



(3-88) 
(3-89) 
(3-90) 

(3-91) 
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Y + X6 + 2X6 






Y + p 



~2 

X + 



(v 



x 2 ~ 
+ p) + Z 



3/2 



= 



(3-92) 



Z +• 



~2 . /~ _\2 ~2] 3 / 2 



X + 



(y +P y 



+ Z' 



(3-93) 



where the dot over the symbol now refers to derivatives with respect to f. It is now 
assumed that the distance between the two vehicles is small compared to the semi- 
major axis of the vehicle at the origin of coordinates. As usual, the nonlinear terms 
are expanded in terms of powers of the coordinates, retaining linear and quadratic 
terms. When this expansion is performed, 



X - Y6 - 2Y<9 + 



1 



X-3XY = 



(3-94) 



y + xe + 2xe 



v^ 



2p 4 



2Y 2 + Z 2 



(3-95) 



z +A - 3^X= 

p 6 p 4 



(3-96) 



For orbits that are nearly circular, the dependence of the angular speed and the 
radius to the center of coordinates upon the independent variable, time, can be taken 
into account by expanding 6 and r s in terms of the eccentricity and the mean 
anomaly. This expansion is done by taking the time derivative of the equation-of-the- 
center (eq. (1-99)). If the equation is truncated at three terms on the right-hand side 



i = M + 2e sin M + | e 2 sin 2M 
4 



(1-99) 



then the ordinary time derivative is 



But 



= M(l + 2e cos M + - e 2 cos 2M + 



a/2 



M = 



S r 
& e e 



(t - t ) 



so 



M 



^ e r e 



a/2 



Hence 



is r \ /r \ / 
f 6 e e 



r 3 

s 



-^■j f 1 + 2e cos M + - e cos 2M + . . . 
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or in nondimensional form 



= ( 1 + 2e cos M + 1 e 2 cos 2M + . . . 



(3-97) 



where, from equation (1-67), 

e 1/2 r 
i3/2 V 1 l o) 



M = 



or from equation (3-89), M is expressed in nondimensional form as, 



M = T - T~ 



Also, 



p = 1 - e cos(t - t ) + 1- j 1 - cos 2(t - t ) 



+ . 



(3-98) 



(3-99) 



Solutions to equations (3-94), (3-95), and (3-96) are to be sought by assuming a 
first-order solution plus a small differential correction for each term. Thus, let 



X — X-j -f- Xo 



Y = Y 1+ Y 2 



Zi — /jA -f - Zn 



(3-100) 
(3-101) 
(3-102) 



If the nonlinear terms are omitted and the origin is assumed to be in a nearly circular 
orbit, equations (3-94), (3-95), and (3-96) revert to the unintegrated form of the 
Clohessy -Wiltshire equations. In nondimensional notation 



X l " 2Y i = 



Y x + 2X, - 3Y, = 

' J l T ~1 



1 " Ji l 
Zi + Z, = 



(3-103) 
(3-104) 
(3-105) 



If, at time t=0, (t = 0),then 



x 1 = 


x 


V 


Yo 


z 1 = 


V 



(3-106) 
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and 



*i = c 



Y l = Y o f 



Z\ - Z 



(3-107) 



The solution^ of equations (3-103), (3-104), and (3-105) become, respectively, 

Xj = 2(x - 3Y )sin t - 2Y Q cos r + 3(2Y - X q )t + (x q + 2Y Q ) (3-108) 

Y t = Y x sin t + 2(5 - 3Y )cos 7 + 2(2Y - X Q ) (3-109) 



5^ = Z< sin t + Z Q cos t 



(3-110) 



It can be seen that these equations are exactly analogous to equations (3-29), (3-30), 
and (3-31). Substituting these solutions back into equations (3-94), (3-95), and (3-96), 
differential equations are obtained for X 2 , Y 2 , and Z 2 . These equations are simpli- 
fied by retaining only the larger quantities such as quadratic terms in the coordinates 
~2~2~2 ««~~ 

X^ ,Y^ ,Zj , etc., and the linear terms in the eccentricity eXj,eYj,eZj_, etc., and 

neglecting the smaller terms. The resulting differential equations are 



X 2 - 2Y 2 



SXjYj + e 



(4Y X + X^cosfr - t ) - 2Yj sinfr - t ) 



Y + 2X 9 - 3Y 9 = JX/ - 2Y, + Z * + e 



£ 2 



2V A 1 



+ 2Xj sinp - t q \ 



(l0Y 1 - 4X 1 )cos(t - t \ 



Z 2 + Z 2 - 3Y 1 Z 1 



3eZ.. cosrr 



>) 



(3-111) 



(3-112) 



(3-113) 



The initial conditions for X 2 ,X 2 , Y 2 ,Y 2 , and Z 2 ,Z 2 are all zero because the gen- 
eral initial conditions have been satisfied by X., Y.,, and Z.. Since the above set 
of equations is a set of linear equations with constant coefficients and known forcing 
functions on the right-hand side, the determination of the solution is straightforward 
but involves considerable effort. The technique is the same as London's with, how- 
ever, still more additional terms. 

Anthony and Sasaki carried out their original analysis by breaking their solution 
into two parts, one dependent on the eccentricity of the orbit and one not dependent on 
the eccentricity. The result is that the part which is independent of the eccentricity is 
formally equivalent to London's solutions in the limiting case of zero eccentricity and 
that the effect of the eccentricity shows up as a second part of each constant. 
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The solutions to equations (3-111), (3-112), and (3-113) are, respectively, 
X 2 = (a ° + eA e ) + (Aj + eA^T + (a 2 ° + eA 2 e Wn t + (a 3 ° + eA 3 e )cos f 
+ (a 4 ° + eA 4 € ^sin 2f + (a 5 ° + eA 5 e )cos 2t + (a 6 ° + eAg € )? sin J 
+ (a 7 ° + eA 7 e )~ cos f (3-114) 

Y 2 = (b ° + eB € ) + (Bj + eBj 6 )? + (b 2 ° + eB 2 e )r 2 + (b 3 ° + eB 3 e ]sin t 
+ (b a ° + eB 4 e )cos f + ^B 5 ° + eB 5 e )sin 2t + (b 6 ° + eBg e )cos 2f 
+ (b 7 ° + eB 7 e lr sin f + (b 8 ° + eBg^r cos r (3-115) 

Z 2 = (c ° + eC e ) + (Cj + eCj^sin t+ (c 2 ° + eC 2 e )cos t + (c 3 ° + eC 3 e )sin 2t 

+ (c 4 ° + eC 4 e )cos 2t+ (c 5 ° + eC 5 e )TSinT+ (c 6 ° + eCg^rcos t (3-116) 

where the coefficients A^° and A^ e , etc., are given by 



A ° = 3(| Y Y - X Y + X Y + I Z Z ) (3-117a) 



o / 9 11 ~ 2 ^2 1^.2 1^,2 i~2 .i ~ ~4.\ 

A ] °=3\X 2 + fY +2X + |y o +|Z +|Z - 13X Y + XoY J (3-117b) 



~ a a ~ a 2 ^2 ~2 ~2~2 ~2 

A 2 ° = -3XoY + 36X Y - 10X o - 2Y - 30Y o - 3X Q - Z - 2Z (3-118c) 



A„° = -6Y Y + 2X Y - 3X Y - 2Z Z (3-118d) 



n ^2,^2 z ~ q^2i~ 2 1~2 
A 4 °=-X + lY + 3X Y -|Y + Iz -iz o (3-118e) 



A 5 ° = X Y - | Y Y + I Z Z Q (3-118f) 
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A 6 ° = -6Y Y + 3X Y 



(3-117g) 



A 7 ° = -21X Y + 6X + 18Y Q 2 



(3-117h) 



o 21 ~ ? 3 ~ 2 3~2 ~ 2 ~~ 3~ 2 3~2 
B = y Y o " 2 Y ° + 2 X ° + 3X ° " 12X « Y o + 4 Z ° + 4 Z ° 



B^ = -6X Q Y + 3X X + 6X Y - 12Y Y 



~ 2 



9 



B 2 =-18Y + 18X Y -|X 



o 1 



2 ^.2 



3* = 



Br 



9,"; 2 .1~ 2 1~2 



- Y + 2X - 6X Y + - Y + - Z - - Z Q 



~ ~ ^,2 „ o 
-21X Y + 6X + 18 Y 



Bg - 6Y Y - 3X Y 



B 3 " = 12Y Y - 7X Q Y - 3X X + 6X Q Y + Z Q Z 



n ~ 2 ~~ ~ 2 ~23~2l~2~2 V 

B 4 ° = 2Y + 18X Y -5X - 15Y Q - § X Q - | Z Q - Z Q ? 



B 5 - 2x o Y o " 3Y o Y o " 2 Z o Z o 



(3-118) 



C o° = o l Y o Z o + 2 *o z o - 3Y o z oJ 



(3-119a) 



Cj - Y Z - XqZ + 3Y Z 



(3-119b) 
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Cg = -2Y Z - 2X Z + 3Y Z 
C 3 ° = - |(y Z + 2X Z - 3Y Z ) 

C 4 = - 4 2X Z - 3Y Z - Y Z j 

C 5 = 2 \ 4Y ° Z ° " 2X o Z cV 
C 6 ° = - -(^4Y Z - 2X Z ^ 

A Q € = ( -3X + I Y Q j sin f + (x - | Y jcos r Q 

A x € = -^3X + 3Y jsin t + (-3X .+ 15Y jcos t 

A 2 e = 3(X + 2Y )sin f Q - 12Y cos t q 

A 3 € = ( 6X - 8Y )sin t + ^-X + 2Y jcos t 

A 4 £ = - § Y sin t + 1(2X0 - 3Y )cos r 

A 5 € = |(-2X + 3Y )sin t - | Y cos r 

A/ = -3(x - 2Y )sin t 
A 7 e = -3|^X - 2Y Jcos f 
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(3-119c) 



(3-119d) 



(3-119e) 



(3-199!) 



(3-119g) 



(3-120) 



B o € = "( 2X o + 3Y )sin r + (-4X + 13Y ) 



COS T„ 



B, =0 



B 2 =° 



Bg = -X sin t q - 2Y cos t 



B 4 £ = 2^X + 2Y )sin t + (2X - 10Y o )cos t q I 



(2X - 3Y )s 



B 5 C = {2X - 3Y Jsin t + Y cos r 



B c 



-Y sin t + (2X - 3Y )cos t q 



B 7 e = (3X - 6Y )cos r 
B 8 e = (-3X + 6Y )sin : ? 



J 



(3-121) 



Co 6 - 



3 ~ ~ 3 ~ 

- Z sin t - - Z Q cos t 



(3-122a) 



Ch 6 = -Z sin t - 2 cos t c 



(3-122b) 



Cg = 2Z sin r + Z cos t 



(3-122c) 



C 3 € = I Z sin ~o + 1 Z cos 7 



(3-122d) 
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C 4 e = - I Z sin f + 1 Z cos T (3-122e) 



C 5 e = (3-122f) 

C 6 e = (3-122g) 

It can be seen from the above two examples of perturbation solutions to relative - 
motion equations that the complexity of the solutions gets out of hand quite rapidly if 
one is to try for high accuracy. The question can then be legitimately asked: what 
advantage can one expect to gain from these solution? The answer, of course, is long- 
term physical insight into the nature of orbits. One often wishes to ask the question: 
will one satellite in the near vicinity of another remain there or drift away and, if so, 
how fast will it drift away? Or, alternatively, where is the best location to place an 
experiment package in the vicinity of, say, a manned space station, so that the longest 
periods of undisturbed experimentation can be carried out? Insight into the answers 
to these questions can be gained by examining the coefficients of the solutions. The 
secular terms will yield the long-term drift rate while the coefficients of the harmonic 
terms are useful for assessing short-term effects. In addition, the determination of 
trajectories which lead to the smallest secular variations give one of the best means 
for establishing initial coordinate locations for experiment packages. It can also be 
seen that orbits which lead to periodic returns to the origin can be studied by this 
means. Such orbits could, of course, be quite useful for recapturing experiment 
packages or, alternatively, for keeping ejected packages away from the origin if so 
desired, thus removing any possibility of future interaction of package and station. 

It can be seen that studies of this nature are highly dependent on the accuracy of 
the secular terms, and a good premium should be placed on getting these terms cor- 
rectly. It will be found in the next chapter that there are advantages along these lines 
to be achieved if rectangular coordinates with their complexity are discarded in favor 
of nonrectangular systems, for just this reason; that is, in some directions, secular 
terms can be simplified and thus lead to better overall results. Rectangular coordi- 
nates are still quite useful in the vicinity of the origin, nevertheless, because of the 
inherent orthogonality of the system. Equations so developed are often simpler to 
implement and visualize, while retaining an acceptable degree of accuracy. It can be 
seen in figures 3-4 and 3-5 for London's equations and in figures 3-6 and 3-7 for the 
Anthony and Sasaki equations that relatively small departure velocities are specified. 
The necessity for the small departure velocities becomes painfully apparent when one 
scales the equations to a small gravitational body such as the moon. In ihe next chap- 
ter an attempt will be made to overcome this difficulty by adopting shell coordinates 
which are themselves curved so as to compensate for the curvature of the orbit. 
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Figure "j-k.- Comparison of London's solution with first-order theory 
and exact theory for a typical earth orbital trajectory in the 
x-direction. 
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Figure 5-5-- Comparison of London's solution with first-order theory 
and exact theory for a typical earth orbital trajectory in the 
y-direction. 
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Figm-e 3-6.- Comparison of Anthony and Sasaki's perturbation solution 
with fir st -order theory and exact theory for a typical earth tra- 
jectory in the x-direction. 
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Figure 3-7.- Comparison of Anthony and Sasaki's perturbation solution with 
first -order theory and exact theory for a typical earth orbital tra- 
jectory in the y-direction. 
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CHAPTER 4 
RELATIVE MOTION IN NONRECTANGULAR COORDINATES 
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L-73-3012 

Apollo Command/Service Module as seen from the Landing Module 
just after rendezvous. The lunar near side is in the background. 
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I CHAPTER 4 

j RELATIVE MOTION IN NONRECTANGULAR COORDINATES 

I 

I 

I In the Previous chapter, relative motion in rectangular coordinates was taken up 

| in some detail. The center of the coordinate frame, however, moves along a curved 

\ path, and the question naturally arises as to what can be gained by the use of nonrec- 

\ tangular coordinate systems. It turns out that, in some cases, such as spherical mov- 

| ing coordinates, nothing in particular can be gained except, perhaps, some modest 

| degree of elegance. If, however, a curved coordinate system centered on the earth is 

I employed, substantial increases in accuracy over long time periods can be attained. 

| This improvement is brought about by various types of more efficient equalization 

j between the dynamic terms and the gravity expansion; that is, curving the coordinate 

| system results in a set of coordinates that is more "natural" to the dynamic physical 

\ situation. This advantage is usually expressed in the final result as a lessening in the 

I number of secular terms necessary to achieve the desired level of accuracy or, alter- 

j natively, as an improvement in the accuracy of the secular term coefficients. This 

j advantage does not mean that on the short time scale curved coordinates are neces- 

j sarily more accurate, although they frequently are, but it does mean that, for large 

\ separation distances or large relative velocities, the curved coordinate systems stay 

j accurate longer. Naturally, the more curved the orbital situation the more one can 

| expect to gain from curved coordinates. This characteristic makes them especially 

\ attractive for smaller gravitational bodies such as the moon. 

| In this chapter, a simple spherical coordinate transformation of the Clohessy- 

[ Wiltshire equations will be taken up first. This transformation provides a simple and 

| useful transition into nonrectangular modes of thought. Next, shell coordinates will be 

I used, taking up, successively, first- and second-order systems. Finally, one of the 

| more successful correction factors which has been tried for the rendezvous equations 

i will be studied. This solution is developed by comparing the first-order rendezvous 

I equations with the known exact solution when both rendezvous and target vehicle are in 

I the same circular orbit. In the course of this development, it will be seen why these 

j systems are good when they are good, and why they fail when they fail. The over- 

j whelming dominance and persistence of the concepts of angular momentum and energy 

\ conservation will also be seen more clearly than before. 



Spherical Coordinates 

The equations used to define the relative motion with respect to a rotating coor- 
dinate system have been presented previously as equations (3-9), (3-10), and (3-11). 
In order to effect the change into spherical coordinates, make the following transfor- 
mations (see fig. 4-1): 
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Figure k-1.- Coordinates employed in describing the motions 
of the space station and shuttle. 
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x = R cos i// cos a ^ 
y = -R cos i// sin 8 > 
z = R sin \]/ ^ 

If this transformation is applied to equation (3-8) the Lagrangian which results is 

1 / • 2 J» 9 • 2 9-9 9 •• .• 

L = - m( R + R z ^ + R'Vcos'^ - 2Rr s cos \p sin -9 + 2Rr g i// sin \p sin 3 



(4-1) 



- 2Rf S S cos i> cos 3 + r g 2 J + m6 s 



2- a _2, 



(Rr s - Rr g )cos \p cos 9 - R 2 3 cos ^ 



+ Rr s i^ sin \p cos 3 + Rr g 3 cos i// sin d 



+ im6 s (r 2 cos 2 i// 



- 2Rr g cos \f/ sin 3 + r g J + 
where 



mg e r^ 



(4-2) 



x 2 + (y + r s ) 2 + z 2 



ill/2 



-( 



R 2 + r g 2 ± 2Rr g cos \p sin d 



.1/2 



The equations of motion corresponding to equations (3-9), (3-10), and (3-11) are 

2\ 



2 

R - Ri// 2 - R(b - S ) cos 2 )// + 



M*.'-^]-. 



cos ip sin 3 



- (r g s + 2r g 6 g )cos \(/ cos Q + 



g Q r 2 R 
e e =0 (4-3) 



|rS + 2R(S - S ) - ReJ cos i// - 2R4/(b - 6 s )sin xp 



+ (r s s + 2r s s jsin 3 



'■(*.' -^-)- f . 



cos 



3 = (4-4) 



Rip + 2R^/ + Rid - 6 S j sin \p cos i// + (r s © s + 2r g e s jsin ip cos 



/ 2 g r \ 



sin \p sin 5 = (4-5) 
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In the special case where the station is in a circular orbit r s and 6 = co are con- 
stants, and equations (4-3), (4-4), and (4-5) reduce to the form 

R-R^ - R(b - w) 2 cos 2 !// + U 2 - -£_§_ r cos ty sin a + -£-£ — =0 (4-6) 

\ V / r v 3 

(R3 + 2R(3 - co)Jcos ^ - 2Ri/$ - co)sin ^ + [co 2 - e e )r s cos -9 = (4-7) 

Ri^ + 2Ri// + R(a - co) sin i£> cos i// - I co - e j- |r s sin p sin 3 = (4-8) 

As before, for the case of a station in a circular orbit, the equations of motion in 
terms of the rectangular coordinate system (eqs. (3-13), (3-14), and (3-15)) may be 
linearized and, therefore, solved in closed form by approximating the gravity differ- 
ence between the two vehicles. Thus 

2 2 

ee Jee ^ . 3 f\ (4 . 9) 



2* o t» o \ 

V S 



This approximation represents a rotating parallel gravitational field rather than the 
spherical field of equation (3-22). The application of equation (4-9) yields equa- 
tions (3-26), (3-27), and (3-28) as was done before in chapter 3: 

x - 2coy = (3-26) 

2 

y + 2cox - 3w y = (3-27) 

2 
z + oTz = (3-28) 

These equations have the solutions: 

x = 2(2 ^ - 3y o )sin cot - 2 ^ cos cot + (ey o - 3 ^jcot + 2 § + x Q (3-29) 

(x \ y x 

2 TT - 3 y J C0S wt + 7J sin cot + 4y Q - 2 -^ (3-30) 

z 

z = z Q cos cot + -£■ sin cot (3-31) 

Inverting these equations to form guidance equations yields the, by now well- 
established, set of equations to rendezvous in time Tj 

x c x (sin co Tl ) + yjsco^ sin co^ - 14(l - cos co^j] 

— .= — — a*. (d-SA) 

w 3coTj sin cot., - 8 (1 - cos cot.) 
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_£ 

W 



2x o( 



1 - COS WT. 



) + y (' 



4 sin ut. 



Zwt* cos WT^ 



3wt.. sin air. 



'1 - cos WT 



l) 



(3-33) 



w tan wt. 



(3-34) 



For simplicity, make the substitution 

A = 3o>T. sin wt. - 8(1 - cos ut-\ 



(4-10) 



In terms of the spherical coordinates R, i^, and -9, equations (3-32), (3-33), and 
(3-34) become, respectively, 

R. / 9 \ O 9 

-^- = f K 2 + K 3 sin d - K 4 sin -9 Q cos S Jcos \p Q + K 1 sin \p Q (4-11) 



*■ 



-£P = (k 1 - K 2 - Kg sin 2 3 + K 4 sin 9 Q cos a Jsin \p Q cos i^ Q (4-12) 



'r 2 

— = Kg sin a cos 9 Q + K^ sin d Q - Kg 



(4-13) 



where 



K, 



Kn 



-COS WT. 



1 



— sin wt, 
A 1 



O 

Ko = —(sin <jot« 



K„ 



K E 



AV 



WT, COS WTj\ 



2 + 2 cos wt^\ 



A* 



COS WT, 



(4-14) 



These then are the Clohessy -Wiltshire equations in spherical-coordinate form. They 
are useful where spherical symmetry about the origin is inherent to the nature of the 
problem as, for instance, if measurements are to be made by a radar located on the 
station. They are not, however, any more accurate than the rectangular version and, 
hence, have seldom been used either in computer studies where the rectangular coor- 
dinates are just as easy to use, or in actual hardware where coordinate conversion has 
been the rule. They do, however, show in an elegant manner the separation of range 
rate from angular rate. 

Shell Coordinates 

It is now informative to look at the form the relative motion equations take in 
shell coordinates. The coordinate system employed in this development (ref. 10) is 
shown in figure 4-2 where r is the projection of a line connecting the center of the 
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planet to the orbital vehicle upon the plane of the reference vehicle. The symbol z 
is normal to this plane passing through the orbital vehicle; y is measured along r 
from the reference vehicle altitude to the projection of the maneuvering vehicle with 
the positive direction upward; and x is measured in a curved arc backward along the 
flight path of the reference vehicle in the plane of the reference vehicle orbit to r. 
The coordinate system rotates about the origin with angular velocity 6. It should be 
noted that this is a left-handed coordinate system as opposed to all of the previous 
coordinate systems which were right-handed. 




>j. Ejection 
\ angle 

i — , 



-^~ 



Direction of 
orbital motion 



Reference vehicle 




■V 



r ^— Orbital vehicle 



Figure It— 2.- Coordinates employed in describing the motions of the vehicles. 

Two assumptions are made about the physical nature of the problem. They are: 
(1) the attracting planetary mass is a gravitational sphere, and (2) the body upon which 
the coordinate system is centered is in a circular orbit. Small departures from these 
assumptions are not considered serious. 
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In a cylindrical coordinate system centered on the planetary body the Lagrangian 



is 



L = i m(r 2 + r 2 2 + z 2 ) - mr 2 wfl + i mto 2 r 2 e e 

2 2 ,/2 2 

Vr + z^ 

Equations (4-1) are converted to shell coordinates by means oi the following 
substitutions: 

y + r s = r 

x = T a e 



Hence 



z = z 



y = r 



x = r c 



z = z 



Then, in shell coordinates, the Lagrangian becomes 



(4-15) 



(4-16) 



(4-17) 



L = ■=• m 
2 



y 2 + z 2 + (y + r s ) 2 (i--co) 



2g r 



e e 



i(y + r s) 



2 2 

+ z 



The differential equations of motion which follow from this Lagrangian are 

x +tl + 221 _ 2wy = 



(4-18) 



(4-19) 



y-(y^m-) 2 - g -7f- 



i + l-f + ' zx2 



-3/2 



= 



g r * 
z + . z 



^) 2 *£ 2 



-3/2 



= 



(4-20) 



(4-21) 



These are the exact differential equations of motion as seen from the orbiting vehicle. 
It can be seen immediately that equation (4-19) is cyclic in the x-coordinate; thus, a 
first integral of the equation of motion in the x-direction is found immediately. The 
integral merely expresses the law of conservation of angular momentum. 



If use is made of the exact orbital expression, 
2 



g o r 



CO 



e e 



(3-16) 
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and equation (4-19) is integrated, the equations of motion become 

' \2 ' ' 



o>r< 



11 + 



k)=« 



(4-22) 



(y + r s)< 



w 



w 



I+ t 



■3/2 



(4-23) 



2 

z + w z 



For convenience, let 



I s 



Y= y_ 
r s } 

z =f 
r s 

T = Wt 



1+^- + ^ 



■3/2 



= 



(4-24) 



(4-25) 



The resulting equations are 
K 



X - 1 = 



(1 + Y)' 



Y - (Y + 1)((X - l) 2 - [(1 + Y) 2 + Z*]~ j 



= 



(4-26) 



(4-27) 



Z + z[(l +Y) 2 + Z 2 



-3/2 



= 



(4-28) 



where the dot over the symbols now refers to the derivative with respect to t. The 
approach taken in obtaining the solutions is the following: Equation (4-27) contains 
terms in X but not in X. The term Z 2 is assumed to be small in relation to 
(1 + Y) 2 and can hence be neglected. For instance, typical dimensional values of z 
and y are on the order of 100 km or less, making Z or Y on the order of 

-i^P- »' 0.107 as compared to 1. Thus, it is possible to cast equation (4-27) as an 
1936 I 

equation in Y and Y by direct substitution of the right-hand side of equation (4-26). 
Equation (4-27) can then be solved approximately for Y as an explicit function of 
time. The solution so obtained may then be used in the solution of equations (4-26) 
and (4-28). 
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In order to establish a value for K, assume that at time equals zero 

T= "") 



x = x > 



Y = Y 



(4-29) 



It follows that K = (x - l)(l + Y ) 2 . Substituting the right-hand side of equa- 
tion (4-26) into equation (4-27) and neglecting the out-of-plane term Z results in 



Y - 



K* 



(1 + Y) 6 (1 + Y)< 



= 



(4-30) 



If the forms of 



(1 + Y) n 



where n = 2 or 3 are expanded, and the terms of the second order and lower are 
retained, equations (4-26), (4-27), and (4-28) become 



X + 2KY - (K + 1) = 3KY 2 



Y + a 2 Y - B 2 = -XY 2 



(4-31) 
(4-32) 



Z + Z = 3YZ 



where 



K=(x o -l)(l + Y o y 



(4-33) 
(4-34) 



a 2 = 3K 2 - 2^ 



K" - 1 



X = -6K* + 3 



(4-35) 



First -Order Solutions 

The first -order solutions to equations (4-31), (4-32), and (4-33) are obtained by 
dropping the second-order terms (setting the terms on the right-hand side of 
eqs. (4-31), (4-32), and (4-33) equal to zero). The Y-equation is solved by inspection. 
Its solution is 



a 2 
Y = a' cos(«t +£')+£_ 



(4-36) 



where primes are used to distinguish first-order -equation integration constants a 
and £ from those of the second-order which will be developed in the next section. 
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The X-equation becomes, upon replacing K + 1 with its equivalent by defini- 
tion X 

X = X - 2KY (4-37) 

Then, substituting the right-hand side of equation (4-36) into equation (4-37) for Y 
and integrating, the X-equation is 

X = (x o - ^V - 2 ^ sin(aT + if) + Mpl sin ?' + X Q (4-38) 

A solution to the Z -equation is 

Z = Z sin t + Z cos t (4-39) 

where Z = Z at r = and Z = Z . 

In order to compute the integration constants from the initial conditions, set 
Y = Y at t = and Y = Y . Then solve for a' and if. This solution can be 
found in a straightforward manner; but, for purposes of making the second-order 
equivalent development which will be undertaken shortly more understandable, it is 
advisable to make the substitution 

u' = a' cos if (4-40) 

Taking the derivative with respect to time of equation (4-36) and applying initial con- 
ditions at t = gives 

Y = -a' a sin ?' (4-41) 

Then by squaring both equations and adding, making use of the identity 
sin if + cos if = 1, a quadratic in a' is obtained 

a 2 u' 2 +Y 2 = a 2 a' 2 (4-42) 

for which the solution 



•4' 2 ^) 2 <«-*» 



where 

n « - v _ §. 

„2 



u'=Y -^ (4-44) 



is obtained by applying terminal conditions to equation (4-36) at r = 0. Hence 



Then it follows easily that 

Jf = cos-^Hl) (4-46) 
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The sign selection on equations (4-45) and (4-46) can be done empirically; for 
instance, a' can be chosen always positive and then the quadrant for if depends 
upon the launch angle 6. Equation (4-36) can be expressed explicitly without the 
epoch angle as 



B 2 Y o 

Y = - £-— cos ar + — 

a 2 a 



sin ar + 



f 



Expanding cos(aT + tf) in equation (4-36) yields 

R 2 
Y = (a' cos if )cos ar - (a' sin if )sin ar + £_ 



(4-47) 



(4-48) 



Comparison of equations (4-47) and (4-48) gives 



Hence, 



a' cos 



a' sin 



tan if = 



_ J3' 



Y r 



-B 2 /a 2 



(4-49) 



(4-50) 



In equation (4-50), since a is a frequency, which is a physical quantity, a nega- 
tive a has no meaning. It can thus be seen that the sign of tan if can be deter- 
mined by using positive values of a . and letting the signs of Y and B 2 determine 
the quadrant. It should be pointed out that B is not necessarily positive due to its 
definition, which permits imaginary values of B. It can be seen, however, that B 
is strictly a function of X and Y and always has the opposite sign to X . The 
reason is that B 2 , assuming reasonably small values of Y Q , is 



x (x - 2) 



It then follows, once a quadrant has been selected for if, that equations (4-49) give 
the sign to be selected for a' in equation (4-45). It is found that a positive selection 
for a' is the correct selection. These characteristics for all possible ejection 
quadrants are given in table 4-1. Thus table 4-1 relates the £' -quadrant and the 
0-quadrant (ejection angle) for the first-order solutions. 

TABLE 4-1.- FIRST-ORDER QUADRANT AND SIGN SELECTION 



x 


Yo 


P 2 


if quadrant 


a' 


<fi quadrant 


+ 


+ 


- 


4 


+ 


1 


- 


+ 


+ 


3 


+ 


2 


- 


- 


+ 


2 


+ 


3 


+ 


- 


- 


1 


+ 


4 
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Second-Order Solutions 
Solution for Y .- Equation (4-32) can be expressed as 

V + a 2 Y + XY 2 = /3 2 (4-51) 
It is advantageous, in order to eliminate the constant )3 on the right-hand side of 
equation (4-51), to transform this equation into the form 

V + p 2 V + XV 2 = (4-52) 

This transformation can be accomplished by a change in variable. Let Y = V + M 
where M is a constant. The transformation results in 

V + (a 2 + 2XM)v + XV 2 = /3 2 - XM 2 - a 2 M (4-53) 

The right-hand side of equation (4-53) thus becomes a function of some still undeter- 
mined constant M. Since the choice of M is not fixed, a value can be chosen which 
makes the right-hand side vanish. Comparing equations (4-51) and (4-52) shows that 

2 - XM 2 - a 2 M = (4-54) 

and 

a 2 + 2XM = p 2 (4-55) 

which can also be written as 



p 2 = \|-15K 4 + 24K 2 - 8 (4-56) 

Solving equation (4-54) for M yields 

2x ^ 4x 2y 

which can also be expressed in terms of p c as 

M =-«i^i (4-58) 

2X 

Equation (4-52) may be recognized as the equation of a one -dimensional anhar- 

monic oscillator. An exact second integral of this equation may be obtained if desired; 

however, as this solution is an elliptic integral, the results are not very useful for 

solving the X- and Z-equations. For this reason an approximate perturbation solution 

is sought. A more complete treatment than that given here using higher-order terms 

can be found in reference 11. The solution V is obtained as the sum of a first-order 

solution Vi and a second-order correction term V^. Let 

V = V x + V 2 (4-59) 

The first-order solution to this equation (by assuming p* » XV) is 

Vj = a cos(pT + (4-60) 
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from which equation (4-59) becomes 
V = a cos(pr + £) + V 2 

Applying these relations to equation (4-52) gives 

2 "2 2 

-p a cos(pr + + V 2 + P a cos(pr + £) + p v 2 



(4-61) 



+ X 



a 2 cos 2 (pr + C) + 2aV 2 cos(pr + ?) + V 2 2 



(4-62) 



(4-63) 



After some cancellation and rearrangement of terms, 

V 2 + p 2 V 2 = -Xa 2 cos 2 (pr + £) - 2XaV 2 cos(pr + £) - XV 2 2 
Omitting terms of higher order than the second (last two terms on right) results in 

(4-64) 



V 2 + p 2 V 2 = -Xa 2 cos 2 (pT + = - i Xa 2 - - Xa 2 cos 2(pr + £) 



Solving the inhomogeneous linear equation in the usual way yields 

2 2 

V 2 = " ^ + ~ cos 2(pr + 
2p^ 6p 2 

Hence, the second-order solution for V is 

2 2 

V = a cos(pr +0 - ^L. + ^_ cos 2(pr + 
2p 2 6p 2 



(4-65) 



(4-66) 



The constants a and £ are integration constants which depend on initial con- 
ditions. This solution is, of course, limited to cases where the first-order assump- 
tion is approximately valid; that is, 



XV « p^ 



or 



X(Y - M) « p 2 



It can be seen from numerical solution that this relationship is usually the case since 
Y seldom exceeds 0.1 for cases which are physically practical and M will be small 
provided the departure velocity is small enough. 

Converting equation (4-66) to Y notation yields 



Y = a cos(pr + £) 

2p2 



1 — cos 2(pr + 
3 F 



« 2 -p 2 
2X 



(4-67) 



The last term on the right is part of a small correction term which takes account 
of the change in energy and hence the change in average altitude between the original 
circular orbit and the new orbit into which the vehicle is launched. 
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Solution for X.- The second-order solution for X is obtained by substituting 
equation (4-67) for Y in equation (4-31) for X and then integrating. The differen- 
tial equation is 



X = X - 2K 



— o 2 9 9 

a cos(pr + + is_ cos 2(pr + ?) _-__ + £_ 

6p 2 2p l iX 2 \ 



+ 3K 



a cos ipr + £) 



+ ^cos 2 2 
36p 4 



4p 4 2p 2 2X 2 x \p 2 */^ x 3p z 



o o o 2 4 

1M_ C0S (pT + - 5LJL C0S( P T + ?) + fi-i COS(pT + Q - ^- COS 2(pT + ?) 



3p<= 



6p" 



2 2 2 

2-2- cos 2(pr + + — cos 2(pT + 

6p 2 6 



(4-68) 



Integrating equation (4-68) with respect to t with A as an arbitrary constant of 
integration yields 

X = A + Bt + C sin(pt + ?) + D sin 2(pT + £) + E sin 3(pr + ?) 



+ F sin 4(pT + £) 



(4-69) 



where 
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A = X - (C sin ? + D sin 2? + E sin 3? + F sin 4?) 



B = K 



2_4 



X t a 2 , p 2 ] 3 ( a 2 . p2 y t a 2 x [ 3g 2 a 2 ^ 19 x 2 a 

K + X + ^2 D 2 o-2 24 o 4 



4X' 



2p Z 



C = -K 



2a ^ 3a(a 2 - p 2 ) , 5Xa 3 
P Xp 2 p 3 



P /\ 4p 2 6p 2 4 p 4 



D=K tfVi-« 2 * ^ 2S 



E = 



F = 



KXa" 
6P 3 



KXV 
96p 5 



(4-70) 



Terms E and F have been found to be negligible in all cases tested but are 
included here for the sake of completeness. The constant A is evaluated on the 
usual assumption that at t = 0, X = X . 

Solution for Z .- The solution for Z is obtained by replacing the right-hand 
side of equation (4-34) by the first -order term for Z and the appropriate second- 
order term for Y to obtain a time -dependent right-hand side. In doing so, of course, 
terms are carried which are higher than the second order of smallness. If only 
second-order terms are retained, however, new end conditions corresponding to a 
and ? would have to be computed as a peripheral calculation. This calculation is 
thought to be unnecessary; therefore, equation (4-33) becomes 

S 



Z + Z = 3(z sin r)/a cos(pr + ?) - 



Xa* 



2p 2 L 



1 - i cos 2(pr + 



si + ei 

2X 2X 



(4-71) 



Let 



Then 



-3Z r 



fXa" 
U> 2 



a 2 -p 2 \_ 



2X 



= K 



Xa 2 Z, 



2(/ 



3aZ = y 



J 



(4-72) 



Z + Z = y cos(pr + ?)sin t + 7] cos 2(pT + £)sin r + k sin t 
The solution of this equation is 



(4-73) 



Z = Z Q sin t + - \yip cos(pr + £) + r;A cos 2(pr + £) + k si 



sin t 



_1 

2PL 



yi^ sin(pr + £) + 2— sin 2(pr + £) + pC 



cos t + Z cos T 



(4-74) 



where 



a 1 2 4 
A = 1 + p + p' 



)//= i +£_ + £_ 
4 16 



C = - 1± sin C - 



^=4 



5p" 



sin 2? 



V 



(4-75) 
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Evaluation of the Integration Constants a and £ 

In order to make effective use of equation (4-67), it is necessary to compute 
values for the two integration constants a and £. It is found from experience that, 
when synchronous or very nearly synchronous orbits are under consideration (that is, 
orbits that return to the same place in space after one complete revolution) , one can 
obtain excellent results by using the values which are obtained from first-order theory. 
For even moderate departures from the synchronous condition, for instance, on the 
order of 5° in launch direction, however, this is not the case, especially for £. It 
will therefore be necessary in most cases to evaluate the integration constants. A 
method for evaluating these constants is outlined in the following. 

Differentiating equation (4-67) and applying the initial conditions yields the 
simultaneous set of transcendental equations 

Y = a cos C - *3_(l - I cos 2?) + ^~ (4-76) 

Xa 2 

Y = -pa sin I - C2- sin 2? (4-77) 

where a and £ are the two unknowns. These two equations can be solved by 
eliminating £ in favor of the two variables u and a where 

u = a cos £ 

After substitution and some manipulation, 

3p 2 3 p : 



^♦"H-St-*" <«-«> 



and 



Solving this set by eliminating a 2 gives a quartic equation in u 



(4-79) 




2.(5p!_2o!. 2Yo ) u 



• 2 2 

| 3 p 2 {a 2 - p 2 ) | Y o | 3 P Y o 

4 x 2 p 2 2X 



(4-80) 



In equation (4-80) it can be seen that u is a function of both Y and X as 
well as Y (through a, p, and X) instead of just X and Y as was the case 
under first-order theory (eq. (4-42)). As a result, the special case where u = will 
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no longer occur independent of Y . It is seen that, when u = 0, the following relation 
is obtained: 



3 P V - P 2 ) , *o 2 , 3 A> 

4 A 2 n 2 + ^T~ 



(4-81) 



or 



° V 4 M >2 2X 



(4-82) 



This relationship is shown in the following sketch for the special case where Y = 0. 
It can be seen that under first-order theory this curve degenerates into a straight 
vertical line. 



u = 



u < 




u < 



u = 



This curve forms the boundary for the selection of the quadrant for £ since 
the sign on u along with the positive selection for a determines the sign of cos ?. 
Since the second-order term in equation (4-67) is a function of cos 2(pr + ?), the 
quadrant is not ambiguous and all four quadrants must be used. Table 4-II summa- 
rizes the second-order sign selection. 

TABLE 4-II.- SECOND -ORDER QUADRANT AND SIGN SELECTION 



Ejection quadrant 


Y 


u 


Sign 
on a 


Quadrant for f 


S Sir 
ir S $ S 2ir 

n S § 2tt 


>0 
>0 
<0 
<0 


>0 
<0 
<0 
>0 


+ 

+ 
+ 
+ 


Fourth 
Third 
Second 
First 
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It can be seen that an equivalent set of circumstances corresponding to the sketch 
and table 4-H exists for each initial value of Y . Thus, in order to map out all the 
possibilities, a series of curves corresponding to u = for each value of Y must 
be computed. This map is done in figure 4-3, for a series of realistic nondimensional 
values. 




Figure k~$.~ Solutions to equation ( lf-82 ) for a series of realistic 



In terms of ejection angle <p, equation (4-82) becomes 



3p% 



taneM^*!!- ^ 2 --^ -^ "° 



x y "* x"x, 



2,v 2 

o 



Ci?*~2\.r 



(4-83) 
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Limitations on the Approximate Solutions 

It will be recalled that one of the conditions implicit in solving the equation of an 
anharmonic oscillator was that p » XV where V was defined as V = Y - M- 



hence, in order for the solution to be valid, 
p 2 » X(Y - M) 



(4-84) 



where Y is inherently nondimensionalized in such a way that, for all practical pur- 
poses, the 'term Y in the expression can be neglected. Thus, the criterion which 
sets a limit on the validity of the solutions in actual practice is 



p » -XM 



(4-85) 



For synchronous speed at ejection angles between 110° and 250°, equation (4-64) 
is not satisfied, as can be seen in figure 4-4. Hence, one would not reasonably expect 
equation (4-67) to describe the physical situation. That this is actually the case can 
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147.0 Synchronous speed 
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Figure h-k.- Comparison of the magnitudes of p 2 and 

and synchronous speeds. 



-AM for Hohmann 



be seen by looking ahead to figure 4-11 for ejection angles of 110° and -110° (250°). 
It is observed that for these cases the error is considerable. At ejection angles of 
slightly larger magnitude than this value, p 2 becomes imaginary (p 2 = ^4X|3 2 + a*). 
Hence, a value for p is undefined. Inspection of equation (4-67) shows that p is 
the frequency term; therefore, no solutions would be expected to exist at all with p 2 
imaginary. This situation occurs for those cases in figure 4-4 for which no solutions 
are shown. For the Hohmann case (that is, apogee at the altitude of the origin with 
perigee at some lower altitude) equation (4-67) is satisfied for all ejection angles and 
p remains real; hence, breakdown does not occur (fig. 4-4). 
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Total -Energy Equation 
Equation (4-30) can be integrated exactly. Multiply through by Y 



YY 



9 • 
KTY 



(1+Y) 3 (1+Y) 2 



Then the integral is 

iL _ 

2 



K' 



1 + C 



x2 1+Y 



(4-86) 



(4-87) 



2(1 + Y)' 
where C is an integration constant. 

Rearranging slightly, C can be identified as the nondimensionalized total 
energy (neglecting the out-of-plane term, of course). Hence, replacing K 2 by its 
equivalent in terms of X and Y 



c-i 



^, ( ltY o) 4 ( i o- 1 ) 2 2 



(1 + YV 



1 +Y 



(4-88) 



An alternative way to derive equation (4-87) is to change the sign on the potential term 
in the Lagrangian, equation (4-18), and then substitute the right-hand side of equa- 
tion (4-22) in this Lagrangian. 



Apogee and Perigee Prediction 
If, at t = 0, Y = Y and Y = Y , 

2 /_._ .\2 2 



-I 



• 2 /• \2 2 

Y o + [ X o ~ 1 ) ~YT 



Y r 



(4-89) 



and remains constant unless disturbed externally. Then equation (4-87) becomes for 
subsequent times 



K' 



(1+Y) 



2 1 +Y 



+ 2C 



(4-87) 



The extremals are then found to be 

-1 ± 



extremal 



^l+2(x -l) 



2 c 



2C 



(4-90) 



The positive value is apogee and the negative value is perigee as measured from 
orbital altitude. 



Solutions in Dimensional Notation 

When equations (4-67), (4-69), and (4-74) are converted back to dimensional 
notation, the following equations result: 
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x = r g [A + But + C sin(po>t + Q + D sin 2 (port + + E sin 3 (port + + F sin 4 (put + g] 

(4-91) 



y = r s ) a cos(pcot + ?) - 



Xa 



1 - - cos 2(pa>t + £) 



2p2L 3 



2X 



(4-92) 



z = r s < t— sin cot + sin cot 



cos cut 



^ cos(pwt + Q+2* cos 2(pcot + + - 
_ 4 4 



2± sin(pcot + + ^ sin 2(pcot + ?) + *^ + £ 



2p " - 4p 

where the constant terms are given by 



2 2 



+ — - cos cot 

1 S 



and 



K = 


x o 

wr s 


1 


a 2 


= 3K 2 


- 2 


/3 2 = 


= K 2 - 


1 



>, 



X = -6KT + 3 



.(4X^ + ^) 1/4 = (. 15K 4 + 24K 2_ 8 )1/^ 



A = j£ - (C sin £ + D sin 2? + E sin 3? + F sin 4?) 



B = K 



2_ + «!L£! + 3fei^ + aV 3« 2 a 2 19 X 2 a 4 " 

XT \ ~ _. *r — — -f- -J- _____ ___™___™ 



cor s K 



4X' 



2p 2 24 p 4 



C= -K 



1* + 3a(a 2 - p 2 ) + 5Xa 3 " 



Xp 



2P-: 



d-k[-- !^^^ 

\ p 4p3 6 p3 4p5y 
3 



E = 



KXa 
6P 5 



KXV 
96p 5 



(4-93) 



(4-94) 



(4 -95a) 
(4 -95b) 

(4-95c) 

(4-95d) 

(4-95e) 

(4-95f) 
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c _- 3 'oA>V».p2^ 

" wr s Up 2 2X , 



77 = 



Xa \ 
2u>r s p2 

3az„ 



air. 



A = 1 + p 2 + p 4 



1// = 1 4. — + £L_ 
^ 4 16 



c = - l± sin ? - p. sin 2? 
P 2p 



(4-95g) 

(4-95h) 

(4-951) 

(4-95J) 
(4 -95k) 

(4-951) 



Equation (4-80) becomes, in dimensional notation, 

2X 2 u 4 /3 | P 2 -q2 2Xy\ L p 2 _ 2g2 2y, 



9P' 



4 \ 2 ' 3p 2 3p 2 rJ 



2X 



,3 pV- P 2 ) , K 2 ,^0 



4 X2 

With equation (4-79) becoming 



p 2 u> 2 r s 2 2Xr s 



(4-96) 



3p 2 y r 



wr g (3p 2 + 2Xu) 



+ u 



and 



C = cos" 1 H 
a 



Equation (4-83) becomes 



tan 



■ = ± 



3 P 4 ( a2 - P 2 )o' 2 r s 2 3p 4 y w 2 r £ 



X x 



2\x 2 



(4-97) 



(4-98) 



(4-99) 



The parameters p 2 , a 2 , and X 2 are functions of x Q only; thus the quadrant 
of £ changes depending upon launch angle 0, but no longer simply at ir/2, vr, 
37r/2, and 2it as was true for the first-order solutions. Figure 4-5 is a plot of x 
as a function of for u = in dimensional terms which apply to the moon. The 
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Figure I4--5.- Plot of u = as a function of x and ejection angle ( 



figure can be used, therefore, to determine for the known value of x and 



the 



proper quadrant for £. As in first -order theory, a is always selected positive. 

Physical Interpretation of Terms 

The physical meaning of several of the terms in these solutions can be seen if 
the first-order solutions and second-order solutions are compared. 

y -equation .- The most important equation from a theoretical standpoint is the 
y-equation. It will be best to begin by comparing equation (4-36) (first-order solu- 
tion) with equation (4-67) (second-order solution). Thus 



Y = a' cos(aT +£')+£_ 
a 2 



(4-36) 



as compared to 



Y = a cos(pr + £) 



Xa 
2p 5 



1 - i cos 2(pr + 



a 



2X 



(4-67) 



In first-order theory, a sinusoidal term is added to a constant which is a function of 
initial velocity in the x-direction. Hence, the constant merely reflects a change in 
energy produced by ejection from the reference orbit and is a function of the speeding 
up or slowing down of the orbital vehicle's angular rate. Superimposed upon this 
motion is an oscillatory motion above and below the mean altitude, the maximum 
amplitudes of which are apogee and perigee, In second-order theory, the same funda- 
mental characteristics appear. In this case, however, they are modified slightly to 
take better account of the inverse -square nature of the planetary gravity field. The 
form is essentially the same as first-order theory presents, but the constant term is 
slightly modified to become the coefficient of the second term combined with the last 
term on the right in equation (4-67). In like manner, the oscillatory term is present 
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although the coefficient is changed slightly. Hence, these two terms represent the 
same type of motion as was present in the first-order theory. Superimposed upon this 
term is an oscillatory term of twice the frequency of the primary term acting very 
much as a second harmonic. The primary effect of this tc*m is to increase the apogee 
and to decrease the perigee. The function of this term is then to take account of the 
weakening of the. restoring force as the gravity field decreases at longer distances 
from the center of the planet. 

Another way of viewing this phenomenon can be seen in figure 4-6 which shows 
the acceleration of the lander relative to the reference orbit in terms of exact first- 
order and second-order theory for a typical case. It can be seen that acceleration is 



Second-order solution 
Exact solution 
First-order solution 



y, m/sec 




a 1 = 0.7273 y Q = 128 m/sec 

P 2 = 0. 0908 u> = , 82 x 10" 3 r ad/sec 



-320 -240 -160 -80 80 160 240 320 x L0° 

y, m 

Figure 4-6.- Acceleration of the orbital vehicle as a function of the 
departure distance from the reference orbit for a typical lunar 
orbit . 

a linear function of displacement over a reasonable altitude range under first-order 
theory. Under second-order theory, the acceleration is quadratic and forms a much 
better approximation to the exact conditions. In figure 4-6, p 2 /a 2 constitutes a 
change in mean equilibrium altitude and shows that, in general, the point of equilibrium 
is displaced, in this case downward. It is pointed out that this change will occur even 
for orbits of the same total energy as the reference orbit since orbits of the same 
period but different eccentricity will have a different mean altitude. Under second- 
order theory this mean displacement is defined a little more precisely by 



Xa 



2X 



lp l 



(not shown in fig. 4-6). It can be seen that, under the special condition that X = 0, 
this relation reduces to fr/a 2 as in first-order theory. 
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Apollo 15 stibsatellite shortly after launch from the Command and Service Module. 
This is the first such vehicle to be launched from a spacecraft in lunar 
orbit . 
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x-equation.- Comparison of the x-equation, equation (4-69), with first-order 
theory, equation (4-38), shows the first-order equation contains three main terms: a 
constant, a secular term, and a sinusoidal term. From the physics of the situation, it 
is clear that the phenomenon which was manifested as a constant mean separation in 
the y-equation, equation (4-67), is a time -dependent linear mean-drift rate in the 
x-equation, equation (4-69). In general, the ejected vehicle would be expected to drift 
away over an interval of several orbital periods due to its difference in energy and 
angular momentum. For a synchronous orbit, cf course, the coefficient of the secular 
term would be zero; thus, this term controls the rate of drift. 

In other respects the form of this equation is about the same as is obtained for 
the y-equation, that is, a series of higher-order trigonometric terms. 

z -equation.- In the z -equation, first-order theory predicts simple harmonic 
motion; so again there is a linear restoring force. This theory then assumes that 
altitude y does not have any appreciable effect upon the out -of -plane mode z. 
Although in trajectories of practical interest, it is quite true that z does not have 
any marked effect upon y because the orbits are r^rly coplanar, the converse of 
this statement is not true; and it is found that, for o^uits of moderately large eccen- 
tricity (for instance, on the order of 0.1), y-coupling into the z -equation can be very 
significant relative to the total, but admittedly small, amplitude in z. In fact, for the 
typical synchronous lunar orbit used for the numerical calculations in this book, the 
coupled term in y and z amounts to one -third of the value of the pure z-term. 
Consequently, the second-order terms in equation (4-71) represent coupling terms 
with y. 

Test Cases 

In order to test the properties of the solutions, the exact differential equations, 
first-order solution equations, and second-order solution equations have been pro- 
gramed on a digital computer. The test orbital conditions were taken as those for 
ejection from a 200 -kilometer circular orbit about the moon. The period of the orbit 
is 7600 seconds. Two particular trajectories are shown in detail. These trajectories 
were a synchronous orbit and a Hohmann case. Both trajectories had a perigee point 
20 kilometers above the surface of the moon. Such orbits are suitable either for 
reconnaissance or for landing, and both orbits were of interest for the Apollo missions. 

Synchronous case .- The chief advantage of selecting the synchronous orbit is 
that, since it returns to its initial relative position after one revolution, the investiga- 
tor is in a position to interpret errors in the trajectories in terms of the particular 
terms in the solutions to which they are due. Such a trajectory is shown as a time 
history in figure 4-7 and in terms of x and y in figure 4-8. It can be seen that the 
second-order solutions predict the time of occurrence of periapsis much fetter than do 
the first-order solutions. 

Ejection is in an upward direction so that apogee is reached at approximately 
one -quarter orbit and perigee at three-quarters orbit. Time histories showing the 
main results of this study are shown in figure 4-7. Particular points of interest are 
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figure 4-7.- Time history for a 200-kilometer synchronous lunar orbit showing 
exact, first-order, and second-order solutions. 
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Figure k-8.- Variation of x with y for a 200-kilometer synchronous lunar orbit 
showing : "the- respective positions of the orbital vehicle as seen from a vehicle 
in reference orbit. 

(1) The error at perigee in these cases for the second-order solution is 5 kilom- 
eters. The second-order solution is, therefore, an improvement over the first-order 
solution where the error is 20 kilometers and seems to indicate that the equations are 
indeed useful as a prediction method for lunar landing or reconnaissance vehicles. 

(2) The second-order y-equation is in error timewise at the end of one orbit by 
200 seconds. Since the period for this orbit is about 6800 seconds, this error is con- 
sidered to be relatively large. The x-equation is in error timewise by the same 
amount. The reason for this time error is found in the nature of the anharmonic 
oscillator equation employed in solving for y. It can be shown that, when the equa- 
tion is restricted to second-order terms, the period is not a function of amplitude. If 
terms of higher order had been carried, the period would be a function of amplitude 
and the time error would be considerably smaller. 

Figure 4 -8 is a plot of the variation of x with y for the synchronous orbit of 
figure 4-7. Hence, this figure shows the position of the ejected vehicle as it would be 
seen from the ejecting vehicle if the out-of -plane motion is disregarded. It is to be 
observed that, in spite of the time error in both x and y, the spatial agreement 
between the exact and approximate second-order equations is very good. 

Hohmann case .- Figures 4-9 and 4-10 show essentially the same information as 
figures 4-7 and 4-8, but for the Hohmann case. First -order data have not been 
included in these cases. The agreement between exact and second -order solutions is 
better for the y -motion than it is for the synchronous orbit since the ejection velocity 
here is much smaller. The data were carried for two complete orbits in order to show 
more clearly the nature of the error buildup which occurs at the end of thyp first orbit 
in the x-equation. It is observed that the approximate solution is better over certain 
portions of the trajectory than over others and that the error is most significant toward 
the end of each orbital period. This result was also obtained in chapter 3 for a rectan- 
gular coordinate system and suggests that the cause of the error is something which is 
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Figure 4-9.- Time history for a 200-kilometer Hohmann transfer to near the 
lunar surface showing exact and second-order solutions. 
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Figure 4-10.- Variation of x with y for a Hohmann transfer from 
200-kilometer trajectory showing successive positions of the 
vehicle as seen by an observer in the orbiting vehicle in the 
reference orbit. 

retained by both approaches. Figure 4-10 shows a gradual departure from the exact 
solutions for each successive orbit because of the secular term of equation (4-69). 

Other launch angles .- In figure 4-11 a series of trajectories is shown for differ- 
ent ejection angles <t> spaced at 20° intervals. The ejection speed is the same as for 
a synchronous orbit. None of these are synchronous orbits, however, because of the 
direction in which ejection takes place. For the sake of clarity, only the exact and 
second-order solutions are shown. 

Figure 4-12 contains the same information as figure 4-9, but the speed is the 
same as that for a Hohmann transfer. Because of the lower ejection velocity, these 
trajectories are in much better agreement than those of figure 4-11. 

In figure 4-11 no approximate solutions are presented over a range of ejection 
angles from -120° to +120°. In addition, the solutions which are presented for ±110° 
are considerably in error. A fundamental assumption in solving the anharmonic 
oscillator equation is that the first-order term is large in comparison with the second- 
order term. For the relatively high ejection speed considered for these cases, the 
assumption is violated over this range of ejection angles. A similar situation is not 
encountered for the lower launch speed of figure 4-12. Mathematically, tne reason for 
this difficulty is that for these trajectories p 2 becomes progressively small in rela- 
tion to the second-order term and then imaginary. When p 2 is nearly as small as 
the second-order term, the solutions are in error and, when p 2 becomes imaginary, 
p is undefined and no solutions exist at all. 
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Figure 4-11.- A series of trajectories for different ejection angles 0. 
Ejection speed is the same as for a synchronous transfer to the lunar 
surface from 200 kilometers altitude. 
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Figure 4-12.- A series of trajectories for different ejection angles 0, 
Ejection speed is the same as for a Hohmann transfer. 

Shell -Coordinate Rendezvous Problem 

It can easily be seen that the solutions to the relative -motion equations expressed 
as equations (4-67), (4-68), and (4-74) are so complicated that reversion to a rendez- 
vous set is difficult if not impossible. It can also be seen that, even if such a rever- 
sion were analytically available, it would most likely not be useful since at this stage 
of development the algebraic complexity begins to outweigh any analytical advantages in 
terms of insight that can be expected. There is, however, a method first suggested by 
Sparrow and Price (ref. 12) which, while not complete in every detail, is # nevertheless 
highly accurate and useful. The method is based upon a comparison of the motions 
predicted by the Clohessy -Wiltshire equations with the motions to be expected if both 
target and rendezvous vehicles are in circular orbits. Since these orbits can be defined 
exactly, it seems logical to compare the exact and approximate results and hence to 
determine just what part of the equations is necessary to define circular trajectories. 



116 




L-73-3014 



Apollo 16 ascent stage above the moon as it approaches the 
Command and Service Module during rendezvous. 
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In this manner a correction factor can be determined which, when added to the conven- 
tional Clohessy-Whitshire set, gives a broad range of accurate applicability. In order 
to keep this problem as simple as possible, it is advantageous to follow the original 
paper and keep the development planar from the beginning. Hence, in what follows, 
the z-coordinate will be neglected with respect to the x- and y-coordinates. 

If equations (4-19) and (4-20) are expanded in the gravity field and only terms of 
the first order are retained, then the following equations result: 

mx - 2m wy = (4-100) 

2 
my + 2mwx - 3ma) y = (4-101) 

It can be seen that these are exactly the same equations as (3-26) and (3-27), except 
that in this case the coordinate system is a shell coordinate system rather than a rec- 
tangular one. Hence, the solutions are the same, namely 




COS tot + I 
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cut 



+ 2 z£- s + r s + ^ < 4 - 102 > 



and 



^ = K _37o\ cos 2c_ 4Yo_2xo 

r s \ wr s r s/ wr s r s wr s y 

where q x and q represent the aggregate of all higher -order terms in the x- and 
y-equations, respectively. Equation (4-22) is, however, an exact relation and, if both 
target and rendezvous vehicles are in circular orbits, y is a constant (the altitude 
difference between the orbits) and equation (4-22) becomes 

A' 1 t 1+ k) =K (4 " 104) 

This equation can also be solved for x yielding the exact expression for circular 
orbits. If at t = 0, x = x and y Q is a constant, K is given by 

v 2 

(4-105) 



(4-106) 




Then equation (4-104) reduces by cancellation to 

x x o 

wr s wr s 

With w canceled this equation integrates directly to 

x x o 

— = — t + Constant 

r s r s 
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or if, at t = 0, x = x , the constant can be determined and 



: _ x o x o 

r r 

s s s 



(4-107) 



The physical explanation for this equation is quite simple. If the difference in 
orbital altitude is held constant, the separation of the vehicles depends solely on the 
difference in their angular rates. This is precisely what x /r s is physically; so the 
separation ajigle is strictly a linear function of time and depends on the initial angular 
separation as an epoch angle. A comparison of equations (4-102) and (4-106) shows 
that for the special case of circular orbits 
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(4-108) 



However, since q is a small correction factor in equation (4-102), it will be suffi- 
cient to drop the smaller parts of q x . For nearly circular orbits obviously y Q ~ 
and for small values of wt will practically cancel anyway. Hence, 
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(4-109) 



If in equation (4-23), y and z are set equal to zero and the equation is solved for 
x /wr g , one obtains for a circular orbit with fixed y 

\-3/2 

(4-110) 



wr c 



= 1 




Then substituting this value back into equation (4-109) for x /wr s the correction fac- 
tor for the x-equation is found to be approximately 
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It is convenient to write this equation in the functional form 
q x = 4Q(y Q )(sin wt - wt) 



where 



^°H^--K] 



3/2 



(4-111) 



(4-112) 



(4-113) 



will be referred to as the "Sparrow Q factor." 

In the y-direction the same situation is obtained. If both orbits are circular, 



^- = ~ !a Constant 



r r 

s s 



(4-114) 
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Then a comparison of equation (4-114) with equation (4-103) shows that 



/2x 3y \ y 2x„ 3y 

q = - — 2. ° cos wt - —2- sin art + — °- 2 (4-115) 

y \ wr s r s/ wr s wr s r s 

As before, the term in y is assumed to be negligible and hence vanishes. Then 

q = .2(| ^2 - -f°_)(i - cos wt) (4-116) 

or by use of equations (4-113) and (4-110) 

q y = -2Q(y )(l - cos wt) (4-117) 

With these correction factors the first-order Clohessy-Wiltshire equations can be 
upgraded in the case of planar motion to the form 

2L = 2 fc°- - ^sin wt - 2 I°-cos art + f 6 ^ - ^°W 
r s \»r s r s/ wr s 



v x 
+ 2^ + ^2 + 4Q(y )(sin wt- wt) (4-118) 

s s 



r s \ wr s r s / 



y 4 y 2x o 

cos wt + sin wt + — - 



- 2Q(y )(l - cos wt) (4-119) 

With considerable cancellation in the Q terms, these equations lead to the equivalent 
equations 

_ x x 2y 

2L = _JL w t + _° + — °.(l -cos wt) (4-120) 



JL - _° + _°_ sin wt (4-121) 

r s r s wr s 

If it is desired to achieve interception at some time Tj, set x and y equal to zero 
and solve for the velocity needed to achieve these conditions. When this procedure is 
followed, i c as a function of Tj can be expressed as # 

x„ x„ sin wt., + 2y„/l - cos utA 

_£_ = _-2 1 °^ ±1 (4-122) 

wr s r s WT l sin WT 1 
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and y can be expressed as 



wr s 



r g sin wtj 



(4-123) 



Equations of Sparrow and Price 

An alternative guidance -equation set is obtained if Q(yA is retained in the 
form of equation (4-113). In this case, the cancellation between equations (4-118) and 
(4-120), and that between equations (4-119) and (4-121) does not occur so extensively, 
and one can obtain guidance equations by setting x in equation (4-118) and y in 
equation (4-119) equal to zero and solving for the conditions needed to rendezvous in 
some specified time t v The set which results is analogous to equations (4-122) and 
(4-123) except, however, that in this case the resultant equations can be expressed as 
the Clohessy -Wiltshire equations with an additional correction factor. 

In equation (4-118) set the left-hand side equal to zero and multiply through by 
sin WTl . In equation (4-119) set the left-hand side equal to zero and multiply through 
by 2(1 - cos wTj\. Then add the two equations to get 



"x 

ur 



— sin wt, + 
r s 1 r s 



2 6(. 



Ct>T| sin WT 



j - 14 + 14 cos wt.] 



Q(y ) 

A 



fwTj + A\ 



where A is given by equation (4-10) 



A = 3wt., sin wt. 



-.( 



1 - COS WT 



■0 



(4-124) 



(4-10) 



Then substitute the right-hand side of equation (4-124) in equation (4-119) for x , 
again with the left-hand side set equal to zero, to get °' 



wr c 



-2^(1 

r R \, 



COS WT. 



i) + r*( si 



sin wt 



1 



3wt. 



Q (y ) 

' ^ A ft- - 3wtj sin cot 1 \u>t 1 



(4-125) 



These (equations (4-124) and (4-125)) are the rendezvous equations of Sparrow 
and Price (ref. 12). It is of some importance to understand the basic physical assump- 
tions under which these equations were derived. Particularly, they were derived 
under the assumption that in the Q term y Q is zero. This assumption is not very 
disadvantageous for real orbits where both rendezvous and target vehicles are in nearly 
circular orbits. Except for this assumption, the equations are found to work extremely 
well so long as both bodies are in the same orbit plane. It is also useful to note that 
the Q-equation (equation (4-113)), degenerates to zero near the target in equa- 
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tions (4-124) and (4-125), thus reducing these equations to the linearized Clohessy- 
Wiltshire form. 

These equations have been tested extensively and have been found to be the best 
purely analytical solutions so far published. The assumptions in the derivation, how- 
ever, do include neglecting secular terms. Thus, these equations will eventually 
break down, and they should be used with caution when rendezvous is to be attempted 
over a long period of time including several orbital periods. 
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CHAPTER 5 

DYNAMICS OF CONNECTED -POINT MASSES 

In chapter 1 the motion of a point mass in orbit was discussed. In chapters 3 
and 4 the relative motion of two point masses was considered in some detail. A logi- 
cal extension of this development is to consider two point masses which are held a 
fixed distance from each other, creating a dumbbell structure. Such a dumbbell in 
orbit under the influence of a gravity field will be studied in this chapter. This dumb- 
bell will then be connected at its center of mass to a second dumbbell through a 
torquing gimbal system. It is hoped that in this manner a simple mathematical model 
can be developed which will allow an analytic treatment with a minimum number of 
approximations and which will exhibit the effects of both gravity gradients and gyro- 
scopic moments. This development is intended to provide a base point for more com- 
plicated analyses. 

Since the dimensions of realistic space vehicles are small in comparison with 
separation distances which have been under discussion in chapters 3 and 4, sufficient 
accuracy for this analysis can be attained by using rectangular coordinates. One 
dumbbell will be assumed to be small and spinning rapidly, thus simulating a gyro- 
scope. The other dumbbell is assumed to be moving slowly and thus simulates the 
characteristics of a cylindrical space station. 

The equations of motion will be found from the Lagrangian subject to the appro- 
priate constraints and then simplified where dropping small terms is appropriate. A 
transformation will then be made on the smaller dumbbell transforming this equation 
to a system more natural to classical gyroscope spin by assuming that precession 
and nutation are slow phenomena. 

Moments Acting on a Simple Dumbbell in a Right -Handed, Moving, 

Rectangular Coordinate System 

The Lagrangian for a single dumbbell in a right-handed, moving rectangular 
coordinate system in a circular orbit is 



T 1 (-2 -2 - 2^ 

L = - m\x + y + z ) 



mw 



1 ™„2 , N 2_2 



(y + r s ) - xyj + - mw (y + r s y + x 



mg r 2 

6 e e 



(y + r s) 



x" + y + r„ + z' 



1/2 



Amw2^x2 + y2 + z 2j (5-1) 



It can be seen that this equation is the same as equation (3-8) with r g set 
equal to zero except for the Lagrange multiplier in the last term on the right which 
acts as a constraint on the distance from the center of mass. Except for these 
changes the assumptions are also the same as in equation (3-8), namely, an inverse- 
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square gravity field for a spherical earth with point masses in orbit. In this case it is 
assumed that both point masses are identical and that motion about the common center 
of mass for the two ends of the dumbell is uncoupled from the orbital motion of the 
center of mass. The constraint specifies that each individual mass is some fixed dis- 
tance p from the common center of mass with p 2 held constant. Thus equa- 
tion (5-1) must be satisfied subject to a constraint (ref. 13) expressed by X, the solu- 
tion for which will allow specification of the tension or compression in the member 
connecting the point masses. 

With the semilength of the dumbbell given by 



p 2 = x 2 + y 2 + z 2 
the equations of motion become 

mx - 2myw - mxw 2 + 



(5-2) 



mg e r e ^x 



x 2 + (y + r s ) + z 2 
2 



3/2 



- 2Xmw 2 x = 



my 



9 m S P r e (y + r s) ? „l 

+ 2mxw - m(y + r B )w* + £_S_i 1— — - 2Xmw z y = > (5.3) 

k 2 + (y + r s ) 2 + z 2 ! 



mz +■ 



mggr^z 



x" + 



( y + r s) 



+ z* 



3/2 



- 2Xmcu 2 z = 



In order to nondimensionalize the foregoing equations for convenience in han- 
dling, divide through by w 2 r s which yields 



mg r 



mx _ 2my _ mx B e 



2X^=0 



w 2 r g wr s r s 




- 2X 



my _ 



(tf-KJ-G 



> (5-4) 



% r e 



mz 



o> 2 r w 2 r 3 
s s 



z/r 



s 



2 / \2 / \2 



— - - 2Xm £- = 
3/2 r s 



J 
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2 g e r e 
Now w = — ■— for one orbit, and so the system can be nondimensionalized in the 
r ° 

s 

normal way for orbits yielding (dropping the unnecessary m) 

-3/2 



X - 2Y - X + X 



X 2 + (1 + Y) 2 + Z 2 



- 2XX = 



-3/2 



2XY = 



Y + 2X - (1 + Y) + (1 + Y)(x 2 + (1 + Y) 2 + z2_ 

Z + Z [_X 2 + (1 + Y) 2 + Z 2 _ - 2XZ = 

with p 2 in equation (5-2) nondimensionalized to the form 
P 2 = X 2 + Y 2 + Z 2 = Constant 
Transposing so as to have only the acceleration terms on the left 

X=2Y+X -X[X 2 + (1 + Y) 2 + Z 2 J +2XX 

Y = -2X + (1 + Y) - (1 + Y)[x 2 + (1 + Y) 2 + Z 2 ]" 3/2 + 

*Z= -z[_X 2 + (1 +Y) 2 + Z 2 ] + 



> (5-5) 



J 



2XY 



2XZ 



The exact moments become 
(1) In pitch 



M, 



pitch = m(XY - YX) 



m T-2XX + X(l + Y) - X(l + Y)[x 2 + (1 + Y) 2 + Z 2 ]" 3/2 
- 2YY - XY + XY[x 2 + (1 + Y) 2 + Z 2 ] " b/2 J 



(2) In roll 

M roll = m(YZ - ZY) 



= m|-Yz[x 2 + (1 + Y) 2 + Z 2 ]" 3/2 + 2XZ - Z(l + 
+ Z(l + Y)[x 2 + (1 + Y) 2 + Z 2 ]" 3/2 j 



Y) 



(5-6) 



(5-7) 



(5-8) 



(5-9) 



126 



(3) In yaw 



M yaw = m(XZ - ZX) 



-3/2 



= mf-Xz[x 2 + (1 + Y) 2 + Z 2 ] - 2ZY - XZ 

_Q /O J 

+ Xz[x 2 + (1 + Y) 2 + Z 2 ] j (5-10) 

It is seen that the Lagrange multiplier cancels out of the moment equations. 
Make a power -series expansion on the gravity term and approximate by dropping 
terms of higher order than the first to get 

|X 2 + (1 + Y) 2 + Z 2 ! ~ 1 - 3Y + Higher-order terms 

Then equations (5-8), (5-9), and (5-10) become 

Mp itch »m[3XY-2<X* + Y*j] (5-11) 

M roll * m(2XZ - 4YZ) (5-12) 

M yaw = m(-XZ - 2ZY) (5-13) 

It can be seen that in the limiting case where the dumbbell velocity is negligible, the 
preceding equations reduce to the usual form of gravity -gradient torques for a rotating 
rectangular coordinate system. The important point here, however, is that the 
dynamic terms are included. These terms allow an adequate treatment of a rapidly 
spinning dumbbell. 

Equations of Motion of a Single Dumbbell System 
Equation (5-6) is 

P 2 = X 2 + Y 2 + Z 2 = Constant (5-6) 

Then the time derivative is 

2PP = 2XX + 2YY + 2ZZ = 
and, taking the second time derivative, 

(XX + YY + ZZ) + (x 2 + Y 2 + Z 2 ) = 
or 

XX + YY + ZZ=-(x 2 +Y 2 + Z 2 ) (5-14) 

This relation can then be used to determine the Lagrange multiplier. Multiply equa- 
tions (5-5) by X, Y, and Z, respectively, and add 
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XX = 2XY + X 2 - X 2 (x 2 + (1 + Y) 2 + Z 2 ]" ' + 2XX 2 



YY = -2XY + Y(l + Y) - Y(l + Y)[x 2 + (1 + Y) 2 + Z 2 ] ' + 2XY 2 
ZZ = -Z 2 [x 2 + (1 + Y)2 + Z 2 ]" 3/2 + 2AZ 2 



) (5-15) 



Then combining these equations by adding the terms containing the Lagrange 
multipliers 



2X 



(x 2 + Y 2 + Z 2 ) = (XX + YY + ZZ) - 2(XY - YX) - [x 2 + Y(l + Y)j 
+ (x 2 + Y 2 + Z 2 )[x 2 + (1 + Y) 2 + Z2]" 3/2 
+ Y [i 2 + (l + Y) 2 + Z 2 ]- 3/2 



When the right-hand side of equation (5-14) is substituted for the first term on the 
right and this equation is solved for X, the result is 

■> _ X 2 + Y 2 + Z 2 XY - YX X 2 + Y(l + Y) 

2(X 2 + Y 2 + Z 2 ) X 2 +Y 2 + Z 2 2(x 2 +Y 2 + Z 2 ) 



+ |[X 2 + (1 + Y) 2 + Z 2 J 3/2 + Y & 2 + (1+Y)2 ^ 2 ] 

2 (x 2 + Y 2 + Z 2 } 



-3/2 



(5-16) 



This expression acts to constrain the length of the dumbbell to a constant value. Thus 
it serves as a measure of the dumbbell tension components in the X, Y, and Z 
directions due to the gravity -gradient and dynamic cross -coupling terms. 

If, now, three direction cosines of the dumbbell orientation are defined by 



cos 6j = 



(X 2 + Y 2 + Z 2 ) 1/2 " P 



and 



Y 

cos 6 = — 

2 P 



cos 6, = — 
3 p 
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It can be seen very readily from equations (5-5) that the nondimensional tension in the 
dumbbell is given by 

T = -XP 

If the gravity term is expanded, X reverts to the approximate form 

■ 2 .2 -2 . 

X + Y + Z XY - YX 1 /.* qv\ 



ro Z2 in 'j ,,2 ji 2 

2[X +Y + Z ) X +Y +Z 



X 2 + 4Y 2 



2(x 2 + Y 2 + Z 2 ) 



(5-17) 



If this value of X is then placed back into equations (5-7) and the gravity terms are 
expanded in these equations also, equations (5-7) become, approximately, 



X - 2Y 1 X 2X ^ - YX) + X (x 2 + Y 2 + Z 2 ) + x(x 2 + 4Y 2 / (5 _ 18) 

X 2 + Y 2 + Z 2 

Y - 2X 1 4Y 2Y ^ ' ^ + y(x 2 + Y 2 + Z 2 ) + y(x 2 + 4Y 2 ; (5 _ 19) 

X 2 + Y 2 + Z 2 

7- 2Z(XY - YX) + z(x 2 + Y 2 + Z 2 ) + z (x 2 + 4Y 2 ) /COA . 

A _ ( 5 _20) 

X 1 + Y 2 + Z 2 

Torque Equations of a Double Dumbbell System 

The equations of motion derived in the foregoing section of this chapter can be 
extended to include the case where the two dumbbells are attached by a suitable gimbal 
system at their centers of mass. To have a perfectly free attachment would be a 
trivial case, since neither dumbbell would be able to sense the other, and the results 
would not be different from the case where each dumbbell is considered separately. 

Consider then, two dumbbells attached by a two-gimbal pivot arrangement (as in 
fig. 5-1) which is free to assume any desired angle but to which torquing motors are 
II attached. The torque applied by either dumbbell becomes a reaction torque on the 

other dumbbell. Let the subscript 1 apply to the large dumbbell which is assumed 
to be moving slowly and which simulates a space station or orbiting satellite, and let 
subscript 2 apply to the smaller, rapidly moving dumbbell which simulates a gyro- 
scope. In order to make the problem as general as possible, it will be assumed that 
the space station can also be controlled by reaction controls and so forth. The equa- 
tions of motion of the space station become, with the summations over all external 
torques , 
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m. 



X 1 Y 1 - YjXi) - 3X lYl + 2 (x^ + YjYj)] = Y T k + J T k (5-21) 



in pitch 



m l[( Y l"Zl - Z lYl) " 2XiZj + 4Y lZ ^] = £ Ti ^ Tl 



in roll, and 



m 1 [(x 1 Z 1 -Z 1 X 1 )+2Z 1 Y 1+ X 1 Z 1 ] = J Tj + £ Tj 



(5-22) 



(5-23) 



in yaw, where the torques produced against the gyroscope are given by Tj i k and 
where external torques (couples) produced by reaction controls, distributed drag, and 
so forth are represented by Tj j k . 



/ 



CONSTANT-SPEED 
GYROSCOPE DRIVE 
MOTOR 




TORQUE 
MOTORS 



TORQUE 
MOTORS 



BASE ATTACHED TO SPACECRAFT-* 



Figure 5-1.- Torquing gimbal assembly. 
The equations of motion of the gyroscope become 




mr 



_X 2 Y 2 - Y 2 X 2 - 3X 2 Y 2 + 2(X 2 X 2 + Y 2 Y 2 )] = 






(5-24) 
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m 2 



(Y 2 Z 2 - Z 2 Y 2 - 2i 2 Z 2 + 4Y 2 Z 2 ) = -£ Ti (5-25) 



and 



m 2 (x 2 Z 2 - Z 2 X 2 + 2Z 2 Y 2 + X 2 Z 2 ) = -^ Tj (5-26) 

If these three equations are combined by elimination of T^ t k > the resultant equa- 
tions are 

X 2 Y 2 - Y 2 X 2 " 3X 2 Y 2 + 2 ( X 2 X 2 + Y 2 Y 2 ) = -(^ 2 )[ X 1 Y 1 " Y 1 X 1 ' 3X 1 Y 1 

+ 2(x 1 X 1+ Y 1 Y 1 )j + ^^T k (5-27) 

Y 2 Z 2 - Z 2 Y 2 - 2X 2 Z 2 + 4Y 2 Z 2 = -(^(y^ - ^ - 2*^ + 4Y lZl ) + ^ T t 



(5-28) 



XpZo - ZoXp HH 2Zpxp *** -"-2 2 



= "(m^A " Z 1 X 1 + 2Z 1 Y 1 +X l Z l) + 4Z T j 



(5-29) 

These equations are general except for the highly accurate gravity expansion and con- 
tain all of the dynamics of interest. 

If it is assumed that there is no thrust from the reaction controls, and it is also 
assumed that the gyroscope is employed in such a way as to maintain the station in an 
essentially motionless condition, the station velocity and acceleration terms are negli- 
gible in relation to the static gravitational torque term. Under these special condi- 
tions the right-hand sides of equations (5-27), (5-28), and (5-29) are dominated by the 
static gravitational torque, and the other terms can be neglected. Other approxima- 
tions are also possible on the left-hand sides of these equations. Small dimensions 
are assumed for the gyroscope; therefore, the static gravitational torque is dropped 
while the velocity and acceleration terms are retained. The resultant simplified set 
of equations of motion becomes 

tH 

11 ..../.. \ /mi\, 

I X 2 Y 2 " Y 2 X 2 + 2 ( X 2 X 2 + Y 2 Y 2 ) = ( 4)( 3X l Y l) < 5 - 30 > 

I Y 2 " Z 2 " Z 2 Y 2 - 2X 2 Z 2 = -(^)(4Y lZl ) (5-31) 

X 2 Z 2 - Z 2 X 2 + 2 Z 2 Y 2 = - ^(X^) (5-32) 
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These equations display two important features. First, gravity-gradient terms 
are still present on the right-hand side of the equations. Second, the gravity -gradient 
terms affect the motion of the gyroscope in a manner that is magnified by the ratio of 
the mass of the station to the mass of the gyroscope. This magnification would seem 
to present a practical limit as to how lightweight space -station control gyroscopes 
should be since the dynamic terms should be capable of balancing the gravity-gradient 
terms for any desired station orientation. 

It is informative at this stage to summarize briefly the assumptions that have 
gone into the derivation of equations (5-30), (5-31), and (5-32). First, it was assumed 
that both station and gyroscope were point-mass dumbbells of fixed length, in orbit, 
connected at their centers of mass but free to move in all directions subject to a 
Lagrangian which implied that the gravity gradient had negligible effect upon the cir- 
cular orbital motion of the center of mass. Next, the two dumbbells were connected 
in such a way that a torque could be applied by one dumbbell to the other. Thus energy 
could be supplied from outside the system by torquing motors while producing no 
change in total angular momentum. 

The equations of motion were simplified by expanding the gravity field as in the 
initial single -dumbbell derivation and subsequently by assuming that the dumbbell 
which symbolized the gyroscope was physically small compared to the station. It was 
further assumed that the gyroscope was employed to maintain the space station in an 
essentially motionless condition which permitted the station velocity and acceleration 
terms to be neglected. This assumption placed a constraint on the lightness of the 
gyroscope relative to the space station. 

Natural Gyroscope Coordinates 

At this stage, both gyroscope and station are expressed in a rectangular coor- 
dinate system. The physical situation is easier to understand, however, if the gyro- 
scopic motion is expressed in a coordinate system which is more natural to the gyro- 
scope. A coordinate system more in agreement with those used in classical mechanics 
is achieved by making the following conversion (see fig. 5-2): 

X 2 = P 2^ cos y cos V + siny cos | sin rj) (5-33) 

Y 2 = P 2 (sin y sin |) (5 . 34) 

Z 2 = P 2^ sin y cos ^ cos ^ " cos y sin V) (5-35) 

where it should be remembered that X 2 , Y 2 , and Z 2 are already nondimensional- 
ized in terms of r so that this set is an exact analogy with equations (1-J-00), y is 
the angle of m 2 in the gyroscope spin plane from a line where this plane intersects 
the horizontal plane, 77 is the angle of intersection of the gyroscope plane with the 
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horizontal plane as measured from the x-axis, and | is the angle between the gyro- 
scope spin vector and the y-axis. This transformation is orthogonal since 



X 2 2 + Y 2 2 + Z 2 2 = P 2 2 



Line of nodes 




>^PIane of motion 
of the dumbbell 



[^Equatorial 
reference plane 



Ascending 
node 



Figure 5-2-- Natural gyroscope coordinates. 
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Under these conditions the pitch equation becomes 

P 2 (sin | cos rj)y + (sin y cos y cos | cos ?? + sin y sin r\)\ + (sin y cos y sin | sin r? 

2 \ •• • 2 / 

- sin y sin ! cos | cos 77)77 - (sin y cos y sin £ cos 77)! + (sin y cos y sin ! cos 77 

+ sin y sin ! cos | sin 77)77 + 2 (cos y cos | cos 77 + sin y cos y sin 77) y| 

+ 2 (sin y sin | cos 77)77! - 2 (sin y cos y sin | cos | cos ?? + sin y sin ! sin 77) y'rj 

( 2 2 2 9 

+ 2 (cos y cos I sin 77 cos 77 + sin y cos y cos | sin 77 + sin y cos y sin ! 

2 2 '\ . / 2 

- sin y cos y cos 77 - sin y cos | sin 77 cos Tjyy + 2 (sin y sin ! cos ! 

2 2 \ • / 2 2 

- sin y cos y sin | sin 77 cos 77 - sin y sin ! cos ! sin 77;! + 2 (sin y cos | sin 77 cos 77 

> 2 2 . . 2 \-l m i / 
+ sin y cos y cos | cos 77 - cos y sin 77 cos 77 - sin y cos y cos ! sin 77 J 77 = i(3X., Y, ^ 



(5-36) 



The roll equation becomes 
™ 2 



(sin I sin 7i)y + (sin y cos y cos | sin 77 - sin 2 y cos 77)*! - (sin 2 y sin ! cos | sin 77 
+ sin y cos y sin £ cos 77)77 - (sin y cos y sin | sin ?7)! 2 + (sin y cos y sin ! sin 77 

- sin y sin | cos £ cos 77)77 + 2 (cos y cos | sin 77 - sin y cos y cos 77) y! 

+ 2(sin y sin | sin 77)77! + 2(sin y sin ! cos 77 - sin y cos y sin | cos | sin 77)^77 
+ 2 (cos y cos I sin 2 77 - sin y cos y cos | sin 77 cos 77 + sin 2 y cos | cos 2 ?; 

- sin y cos y sin 77 cos 77) y + 2 (sin y sin £ cos ! sin 77 cos 7? - sin y cos y sin £ sin 2 ??)! 
+ 2 (sin y cos y cos ! sin 77 cos 77 - sin 2 y cos 2 ! cos 2 ?7 + sin y cos y cos ! sin 77 cos 77 

- COs2ysin2 ^]=-£( 4Y l Z l) (5-37) 
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and the yaw equation becomes 



(cos |)y - (sin y cos y sin |)£ - ^sin y cos £ + cos yjrj - (sin y cos y cos £)| 
- 2 (cos y sin ny| + (sin y sin | cos |jfj^" + 2 (sin y cos y sin |JyT? 
+ 2 (sin y cos y sin | cos | cos 77 - cos y sin | sin r\jy + 2 (sin y cos | cos 77 



- sin y COS y cos 



I sin v)i 



m 2 y X l Z l 



(5-38) 



Owing to the choice of this coordinate system the. time for one complete revolu- 
tion of the gyroscope will be very short in relation to the time required for precession 
or nutation; that is, y is much larger than either £ or 77. Thus, trigonometric 
functions of y can be averaged on the basis that the other parameters are essentially 
constant over any given single cycle yielding considerable simplification in the equa- 
tions of motion. It can be readily surmised that the penalty for this approximation 
consists in ignoring the fast nutation. 

For a single turn of the gyroscope through an angle 2ir, definite integrals can be 

used; that is, 

p2ir p27r 

\ sin x dx = \ cos x dx = (5-39) 



i 



2ir 



sin x cos x dx = 



(5-40) 



I 



2tt 



.. 2 



2ir 



sin x dx = \ cos x dx = 77 
J 



Hence, the average value of the even terms is 



(5-41) 



2 it Jr 



< p2ii is 1 c2t! 

i- sin 2 : , __ 

2rr Jq 2ir J o 2 



cos x dx = 



(5-42) 



and the odd terms vanish. With this spin averaging, the pitch, roll, and yaw equations 
become, respectively, 



1 " 1 1 9 

(sin | cos rj)y + -(sin 77) | - -(sin £ cos £ cos 77)77 + ±(sin t, cos | sin 77)77 



+ (cos I cos r\)yl + (sin | cos 77)77! - (sin | sin rj)'yri + (sin | cos | - 



sin t cos 



- (sin I sin 77 cos 77)77 



m, 



= =r±/'3X 1 Y 



m 



2 V 



1M 



I sin 2 ??)! 



(5-43) 
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1 *• 1 "1 -2 

(sin £ sin rj)y - -(cos 77) | - -(sin £ cos | sin 77)77 - -(sin 5 cos | cos 77)77 

£ Z it 



(cos I sin Tj)y| + (sin % cos 77)777 + (sin | sin 77J77I + (cos |)y 



II 2 2 V" 

+ (sin I cos I sin 77 cos 77)1 - (^cos | cos 77 + sin 77^/77 



m 



1 



) 



(5-44) 



and 



(cos |)y - -(cos I + 1)77 - (sin ij)y| + (sin | cos 4)774 - (sin 4 sin 77)7 



+ (cos I cos 77JI = - i g-(X 1 Z 1 ) 



(5-45) 



The gyroscope is now specified by slowly varying parameters which define the position 
of the spin axis as a function of time. 

Separation of the acceleration terms can be accomplished as follows: Multiply 
the pitch equation by cos 4 cos 77, the roll equation by cos 4 sin 77, and the yaw equa- 
tion by sin |; add the two former and subtract the latter to obtain 



-(sin 1)77 + ri + (sin 77)7 - (cos 4 sin 77)77 

Li 



m l \ 

= jjp^XjYj cos I cos 77 - 4Y 1 Z 1 cos 4 sin 77 + XjZj sin 4 J 



(5-46) 



Then multiply the pitch equation by sin 77 and the roll equation by cos 77 and 
subtract 



- 4 + -(sin I cos 1)77 - (sin 1)777 - (cos 4 cos 77)7 + (cos 77J77 



m 1 

^(SX^sin^Y^cos,) (5-47) 



Finally, multiply the pitch equation by sin 4 cos 77, the roll equation by sin 4 sin 77, 
and the yaw equation by cos 4. Add these to obtain 



7 - (cos 1)77 + (sin 4)4?? + (cos £ cos 77)4 - (sin £ sin 77)77 



m. 



m 



2^ 



(3X..Y., sin 4 cos 77 - 4Y< Z. sin 4 sin 77 - XjZ. cos 4) 



(5-48) 
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These equations represent the motion of the gyroscope when employed to hold a space 
station that is under the influence of a gravity -gradient field stationary. The solution 
of equation (5-46) gives the precession rate, the solution of equation (5-47) gives the 
nutation rate, and the solution of equation (5-48) gives the change in spin speed result- 
ing from conditions on the station specified by X,, Y-,, and Z-,. If other control 
maneuvers are to be made, the more general relationships for station motion and 
acceleration may be required as given in equations (5-27), (5-28), and (5-29). Solu- 
tions of equations (5-36), (5-37), and (5-38) with the more general expressions for 
station effects would give the torque program required by reference back to equa- 
tions (5-24), (5-25), and (5-26). 

Equations (5-46), (5-47), and (5-48) can also be written in terms of moments of 
inertia. If 

h = m 2 P 2 2 
is the nondimensional principal moment of inertia of the gyroscope dumbbell and 

l XY 1 l 1 
h = m l Y l Z l 



h xz = m l X l Z l 
are the respective cross products of inertia of the station, then 
i(sin |)?7 + y| + (sin r\)'y - (cos £ sin 77)77 



31, 41, L 

XY y Z XZ 

— — cos £ cos 77 i — cos I sin 77 + sin | (5-46a) 

h l 2 h 



1-1 -2 / 3 \, 

- £ + -(sin £ cos Dr? - (sin ^)y'rj - (cos £ cos 77)7 + \cos 77J77 
2 & 

3I X 4^ 

= — — sin t? + — — cos 77 (5 -47a) 

l 2 h 



and 



y - (cos £)?7 + (sin Dl'77 + (cos £ cos 77) | - (sin £ sin 77)77 

31 4I i h 

= _2X. sin I cos 77 — sin | sin 77 — cos | (5-48a) 

h h h 

These equations apply either dimensionalized or not since they involve the ratio of the 

inertias as the only dimensional terms. 

137 



Langley Research Center, 

National Aeronautics and Space Administration, 
Hampton, Va., March 19, 1973. 



138 



SYMBOLS 

K 2 
a constant defined as cos y, - — y b 

a constant of integration which describes amplitude of departure from 

reference circular orbit 

,♦- Ko 

b constant defined as — - 

P 

C constant 

Crj drag coefficient 

Cl lift coefficient 

D drag force 

E eccentric anomaly 

g gravitational acceleration 

g earth gravitational acceleration 

H total energy per unit mass 

H s scale height 

h angular momentum per unit mass 

h height above surface of earth 

h a maximum (apogee) altitude attained by vehicle 

hp minimum (perigee) altitude attained by vehicle 

Ij moment of inertia of gyroscope in spin 

Ig = I3 moment of inertia of gyroscope perpendicular to spin axis 

i,j,k orthogonal coordinates 

K integration constant defined in equation (4-34) 
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c nP e s 

K, constant defined as 

1 2m 



C l/e S 

Ko constant defined as 

1 2m 

L Lagrangian, potential energy subtracted from kinetic energy 

L lift force in chapter 2 

L/D lift -drag ratio 

M mean anomaly 

M translation substitution used to solve equation (4-52) and defined 
in equation (4-57) 

m mass of vehicle 

m mass of dumbbell when single or of station -gyros cope combination when 
two dumbbells are present 

P orbital period 

p pressure 

Q(y ) function defined in equation (4-113) 

q dynamic pressure 

R radial distance from space station to vehicle 

r distance from center of planet to orbital vehicle 

r e radius of earth 

r m radius of attracting body 

r s distance from center of attracting body to reference orbit 

S surface area of entry vehicle 

s distance traveled 
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t time 

u dummy variable, a cos £ 

V velocity 

V substitution variable used to solve equation (4-51) 
V c i r circular satellite velocity at surface of earth 

W weight of entry vehicle 

X,Y,Z nondimensional coordinates of shell coordinate system centered on a body 

moving in a circular orbit about an attracting planet or satellite 

x,y,z dimensional coordinates of shell coordinate system centered on a body 

moving in a circular orbit about an attracting planet or satellite 

Y function of a, b, and y defined by equation (2 -50a) 
y function of altitude, e"^ 

a angle of attack 

a,|3,X functions of initial velocity component in x-direction defined 

in equation (4-35) 

/3 decay constant of atmosphere 

T function defined by equations (2-56a), (2-56b), and (2-56c) depending 

on conditions 

y flight -path angle 

Y spin angle of gyroscope dumbbell from line of nodes 

y,r),K constant functions of initial out-of -plane velocity defined in equation (4-72) 

A function of a power expansion in p defined in equation (4-75) 

e eccentricity 

k arbitrary constant of integration which describes epoch angle with 

respect to reference circular orbit 
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J 



jl(y£ 



£(V) function of velocity, -°^[~ - llcos y 

7] angle between line of nodes and x-axis measured in x-z plane 

6 range measured in angle about center of earth 

6 angular coordinate in a cylindrical system with origin at center of 

attracting body and lying in plane -of -reference circular orbit, 
measured from positive x-direction clockwise 

S angle formed by x-axis and projection of R on orbital plane of space 

station 

X Lagrange multiplier 

li gravitational constant 

£ angle between y-axis and angular -momentum vector of gyroscope 

P nondimensional semilength of dumbbell 

p density of atmosphere 

p function of a, j3, and X defined in equation (4-55); hence, a function 

of initial velocity and position in x-direction only 

p nondimensional radius defined by equation (3-88) 

p semilength of dumbbell 

t nondimensionalized or scaled time, wt 

t 1 time to rendezvous 

$ ejection angle 

\p function of a power expansion in p defined in equation (4-75) 

\p angle formed by R and orbital plane of space station * 

0, directional angular velocity 

u) angular velocity (magnitude) 
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Subscripts: 

a apogee or highest trajectory point 

b breakout, or exit, condition 

c computed 

e * earth 

i,j,k orthogonal coordinates 

m moon 

max maximum 

initial condition 

p perigee or lowest trajectory point 

s station position or location of origin 

v vehicle 

1 space station 

2 gyroscope 

Dots over symbols denote differentiation with respect to time or scaled time 
depending upon the section of the book in which they appear. Primes refer to first- 
order-solution quantities. The tilde ~ represents nondimensionalization with 
respect to the semimajor axis. 

For constants, any consistent set of units may be used. In this paper the fol- 
lowing values were chosen: 

g = 9.807 m/sec 2 

r e = 6378.15 km 

r m = 1736.5 km 

r g = r m + 200 = 1936.5 km 

j3 = 0.13666 per km 
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M = 4.8936 X 10 12 m 3 /sec 2 
v cir = 7848.62 m/sec 
p = 1.225 kg/m3 
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